Análisis Comparativo de Redes de Policonsumo

Red de Co-ocurrencia vs Red Bipartita

Autor/a

Amaru Simón Agüero Jiménez

Fecha de publicación

10 de septiembre de 2025

1 REDES DE POLICONSUMO DE SUSTANCIAS EN CHILE

1.1 CONFIGURACIÓN INICIAL

Código
# Función para instalar y cargar paquetes
load_packages <- function(packages) {
  for (pkg in packages) {
    if (!require(pkg, character.only = TRUE, quietly = TRUE)) {
      install.packages(pkg, dependencies = TRUE, repos = "https://cran.r-project.org")
      library(pkg, character.only = TRUE)
    }
  }
}

# Lista de paquetes necesarios
required_packages <- c(
  "tidyverse", "igraph", "tidygraph", "ggraph", 
  "kableExtra", "DT", "plotly", "visNetwork", 
  "corrplot", "RColorBrewer", "patchwork", "dendextend"
)

# Cargar todos los paquetes
load_packages(required_packages)

# Configuración global
options(scipen = 999)
theme_set(theme_minimal())

# ===============================================
# FUNCIÓN DE SIMULACIÓN MEJORADA BASADA EN DATOS REALES CHILENOS
# ===============================================

simular_datos_chile <- function(n_pacientes = 223055, seed = 42) {
  
  set.seed(seed)
  
  # ===============================================
  # DISTRIBUCIONES EXACTAS DE LOS DATOS REALES
  # ===============================================
  
  # Distribución exacta de sustancia principal
  dist_principal <- c(
    "Alcohol" = 0.3507386071,
    "Pasta Base" = 0.3773329448,
    "Cocaína" = 0.1886216404,
    "Marihuana" = 0.0618546995,
    "Sedantes" = 0.0116832171,
    "Opioides" = 0.0044652664,
    "Otros" = 0.0012104638,
    "Inhalables" = 0.0008069759,
    "Hipnóticos" = 0.0006500639,
    "Estimulantes" = 0.0006007487,
    "Anfetaminas" = 0.0005872991,
    "Éxtasis/MDMA" = 0.0003586559,
    "Alucinógenos" = 0.0002645088,
    "Crack" = 0.0002465760,
    "Metanfetaminas y Otros Derivados" = 0.0001613952,
    "Tranquilizantes" = 0.0001613952,
    "Lsd" = 0.0001255296,
    "Heroína" = 0.0000806976,
    "Fenilciclidina" = 0.0000224160,
    "Metadona" = 0.0000224160,
    "Esteroides Anabólicos" = 0.0000044832
  )
  
  # Probabilidad exacta de policonsumo por sustancia principal
  prob_policonsumo <- c(
    "Pasta Base" = 0.8439274767,
    "Alcohol" = 0.4834854411,
    "Cocaína" = 0.8360706391,
    "Marihuana" = 0.7993041966,
    "Sedantes" = 0.6266308519,
    "Opioides" = 0.6817269076,
    "Otros" = 0.6333333333,
    "Inhalables" = 0.8333333333,
    "Hipnóticos" = 0.7103448276,
    "Estimulantes" = 0.6791044776,
    "Anfetaminas" = 0.8396946565,
    "Éxtasis/MDMA" = 0.9875000000,
    "Alucinógenos" = 0.9491525424,
    "Crack" = 0.7454545455,
    "Metanfetaminas y Otros Derivados" = 0.9722222222,
    "Tranquilizantes" = 0.7777777778,
    "Lsd" = 0.9285714286,
    "Heroína" = 0.9444444444,
    "Fenilciclidina" = 0.6000000000,
    "Metadona" = 0.2000000000,
    "Esteroides Anabólicos" = 1.0000000000
  )
  
  # Distribución exacta del número de sustancias
  dist_n_sustancias <- c(
    "1" = 0.2903857793,
    "2" = 0.3405438121,
    "3" = 0.2488579050,
    "4" = 0.1202125036
  )
  
  # Probabilidades condicionales completas
  prob_condicionales <- list(
    "Pasta Base" = list(
      "Alcohol" = 0.6517239741,
      "Marihuana" = 0.5286576527,
      "Cocaína" = 0.3013687237,
      "Sedantes" = 0.0226694865,
      "Otros" = 0.0136397120,
      "Inhalables" = 0.0097664140,
      "Anfetaminas" = 0.0052871706,
      "Opioides" = 0.0022693249,
      "Estimulantes" = 0.0021980372,
      "Crack" = 0.0016633795,
      "Alucinógenos" = 0.0013307036,
      "Lsd" = 0.0011406031,
      "Éxtasis/MDMA" = 0.0008316898,
      "Metanfetaminas y Otros Derivados" = 0.0008316898,
      "Hipnóticos" = 0.0004990139,
      "Heroína" = 0.0004039636,
      "Tranquilizantes" = 0.0003326759,
      "Hongos" = 0.0002613882,
      "Metadona" = 0.0000831690,
      "Fenilciclidina" = 0.0000237626,
      "Esteroides Anabólicos" = 0.0000118813
    ),
    "Alcohol" = list(
      "Marihuana" = 0.2604622031,
      "Cocaína" = 0.2284684408,
      "Pasta Base" = 0.1309788583,
      "Sedantes" = 0.0465015211,
      "Otros" = 0.0214101286,
      "Anfetaminas" = 0.0052406882,
      "Inhalables" = 0.0051000844,
      "Estimulantes" = 0.0045632334,
      "Opioides" = 0.0022880078,
      "Éxtasis/MDMA" = 0.0016872460,
      "Lsd" = 0.0015977708,
      "Alucinógenos" = 0.0011631771,
      "Hipnóticos" = 0.0011631771,
      "Metanfetaminas y Otros Derivados" = 0.0009714446,
      "Tranquilizantes" = 0.0006263261,
      "Hongos" = 0.0005879797,
      "Crack" = 0.0005496332,
      "Heroína" = 0.0003067720,
      "Esteroides Anabólicos" = 0.0001406038,
      "Metadona" = 0.0000511287
    ),
    "Cocaína" = list(
      "Alcohol" = 0.7019941530,
      "Marihuana" = 0.4234069356,
      "Pasta Base" = 0.1497635063,
      "Sedantes" = 0.0355097093,
      "Otros" = 0.0151878877,
      "Anfetaminas" = 0.0101252585,
      "Éxtasis/MDMA" = 0.0068927816,
      "Lsd" = 0.0050626292,
      "Estimulantes" = 0.0040405961,
      "Opioides" = 0.0033275497,
      "Alucinógenos" = 0.0029472583,
      "Inhalables" = 0.0023292848,
      "Metanfetaminas y Otros Derivados" = 0.0018539206,
      "Hipnóticos" = 0.0017350795,
      "Hongos" = 0.0011408742,
      "Tranquilizantes" = 0.0010458013,
      "Crack" = 0.0008794239,
      "Heroína" = 0.0004278278,
      "Esteroides Anabólicos" = 0.0001426093,
      "Metadona" = 0.0000713046,
      "Fenilciclidina" = 0.0000475364
    ),
    "Marihuana" = list(
      "Alcohol" = 0.6352105530,
      "Cocaína" = 0.3405812858,
      "Pasta Base" = 0.2098282235,
      "Sedantes" = 0.0484163224,
      "Otros" = 0.0232659274,
      "Lsd" = 0.0130463144,
      "Éxtasis/MDMA" = 0.0119591215,
      "Inhalables" = 0.0091324201,
      "Anfetaminas" = 0.0081901863,
      "Alucinógenos" = 0.0071029934,
      "Estimulantes" = 0.0068130753,
      "Hongos" = 0.0050735667,
      "Opioides" = 0.0047836486,
      "Metanfetaminas y Otros Derivados" = 0.0015220700,
      "Tranquilizantes" = 0.0009422338,
      "Hipnóticos" = 0.0007972748,
      "Crack" = 0.0007247952,
      "Fenilciclidina" = 0.0002899181,
      "Heroína" = 0.0002174386,
      "Metadona" = 0.0000724795
    ),
    "Sedantes" = list(
      "Alcohol" = 0.4712202609,
      "Marihuana" = 0.2256331543,
      "Cocaína" = 0.1435149655,
      "Pasta Base" = 0.0629316961,
      "Otros" = 0.0376055257,
      "Opioides" = 0.0226400614,
      "Estimulantes" = 0.0072908672,
      "Éxtasis/MDMA" = 0.0072908672,
      "Hipnóticos" = 0.0072908672,
      "Anfetaminas" = 0.0042210284,
      "Inhalables" = 0.0042210284,
      "Metanfetaminas y Otros Derivados" = 0.0023023791,
      "Hongos" = 0.0019186493,
      "Alucinógenos" = 0.0007674597,
      "Esteroides Anabólicos" = 0.0007674597,
      "Metadona" = 0.0007674597,
      "Tranquilizantes" = 0.0007674597,
      "Crack" = 0.0003837299,
      "Lsd" = 0.0003837299
    ),
    "Opioides" = list(
      "Marihuana" = 0.3584337349,
      "Alcohol" = 0.3463855422,
      "Cocaína" = 0.1987951807,
      "Sedantes" = 0.1506024096,
      "Pasta Base" = 0.0793172691,
      "Otros" = 0.0301204819,
      "Estimulantes" = 0.0100401606,
      "Hipnóticos" = 0.0100401606,
      "Anfetaminas" = 0.0090361446,
      "Metadona" = 0.0090361446,
      "Tranquilizantes" = 0.0080321285,
      "Éxtasis/MDMA" = 0.0070281124,
      "Metanfetaminas y Otros Derivados" = 0.0070281124,
      "Hongos" = 0.0060240964,
      "Crack" = 0.0030120482,
      "Fenilciclidina" = 0.0020080321,
      "Lsd" = 0.0020080321,
      "Alucinógenos" = 0.0010040161,
      "Inhalables" = 0.0010040161
    ),
    "Otros" = list(
      "Alcohol" = 0.4111111111,
      "Marihuana" = 0.3074074074,
      "Cocaína" = 0.2148148148,
      "Pasta Base" = 0.1333333333,
      "Sedantes" = 0.0740740741,
      "Opioides" = 0.0370370370,
      "Inhalables" = 0.0296296296,
      "Hipnóticos" = 0.0148148148,
      "Éxtasis/MDMA" = 0.0111111111,
      "Tranquilizantes" = 0.0074074074,
      "Alucinógenos" = 0.0037037037,
      "Estimulantes" = 0.0037037037
    ),
    "Inhalables" = list(
      "Alcohol" = 0.6222222222,
      "Marihuana" = 0.4944444444,
      "Pasta Base" = 0.3222222222,
      "Cocaína" = 0.1000000000,
      "Otros" = 0.0555555556,
      "Sedantes" = 0.0388888889,
      "Alucinógenos" = 0.0111111111,
      "Opioides" = 0.0055555556
    ),
    "Hipnóticos" = list(
      "Alcohol" = 0.4620689655,
      "Cocaína" = 0.2413793103,
      "Marihuana" = 0.2068965517,
      "Sedantes" = 0.1034482759,
      "Pasta Base" = 0.0758620690,
      "Otros" = 0.0551724138,
      "Alucinógenos" = 0.0068965517,
      "Anfetaminas" = 0.0068965517,
      "Estimulantes" = 0.0068965517,
      "Tranquilizantes" = 0.0068965517
    ),
    "Estimulantes" = list(
      "Alcohol" = 0.4104477612,
      "Marihuana" = 0.3582089552,
      "Cocaína" = 0.2089552239,
      "Sedantes" = 0.0895522388,
      "Pasta Base" = 0.0597014925,
      "Éxtasis/MDMA" = 0.0373134328,
      "Anfetaminas" = 0.0223880597,
      "Otros" = 0.0223880597,
      "Lsd" = 0.0149253731,
      "Opioides" = 0.0149253731,
      "Metanfetaminas y Otros Derivados" = 0.0074626866
    ),
    "Anfetaminas" = list(
      "Alcohol" = 0.6335877863,
      "Marihuana" = 0.3664122137,
      "Cocaína" = 0.3129770992,
      "Pasta Base" = 0.1221374046,
      "Sedantes" = 0.0763358779,
      "Estimulantes" = 0.0305343511,
      "Éxtasis/MDMA" = 0.0229007634,
      "Opioides" = 0.0152671756,
      "Alucinógenos" = 0.0076335878,
      "Lsd" = 0.0076335878,
      "Metanfetaminas y Otros Derivados" = 0.0076335878,
      "Tranquilizantes" = 0.0076335878
    ),
    "Éxtasis/MDMA" = list(
      "Marihuana" = 0.6875000000,
      "Alcohol" = 0.5250000000,
      "Cocaína" = 0.4875000000,
      "Sedantes" = 0.2000000000,
      "Opioides" = 0.1125000000,
      "Hongos" = 0.1000000000,
      "Alucinógenos" = 0.0625000000,
      "Estimulantes" = 0.0375000000,
      "Lsd" = 0.0375000000,
      "Otros" = 0.0250000000,
      "Pasta Base" = 0.0250000000,
      "Inhalables" = 0.0125000000,
      "Metanfetaminas y Otros Derivados" = 0.0125000000
    ),
    "Alucinógenos" = list(
      "Marihuana" = 0.6271186441,
      "Alcohol" = 0.5593220339,
      "Cocaína" = 0.3050847458,
      "Éxtasis/MDMA" = 0.1016949153,
      "Sedantes" = 0.1016949153,
      "Pasta Base" = 0.0677966102,
      "Hongos" = 0.0508474576,
      "Estimulantes" = 0.0338983051,
      "Lsd" = 0.0338983051,
      "Anfetaminas" = 0.0169491525,
      "Opioides" = 0.0169491525
    ),
    "Crack" = list(
      "Marihuana" = 0.5454545455,
      "Alcohol" = 0.5272727273,
      "Cocaína" = 0.3090909091,
      "Pasta Base" = 0.1090909091,
      "Anfetaminas" = 0.0181818182
    ),
    "Metanfetaminas y Otros Derivados" = list(
      "Alcohol" = 0.6111111111,
      "Marihuana" = 0.3888888889,
      "Cocaína" = 0.1944444444,
      "Anfetaminas" = 0.0833333333,
      "Estimulantes" = 0.0833333333,
      "Hipnóticos" = 0.0833333333,
      "Pasta Base" = 0.0555555556,
      "Sedantes" = 0.0555555556
    ),
    "Tranquilizantes" = list(
      "Alcohol" = 0.4722222222,
      "Marihuana" = 0.2500000000,
      "Cocaína" = 0.1111111111,
      "Estimulantes" = 0.0833333333,
      "Opioides" = 0.0833333333,
      "Otros" = 0.0555555556,
      "Hipnóticos" = 0.0277777778,
      "Lsd" = 0.0277777778,
      "Sedantes" = 0.0277777778
    ),
    "Lsd" = list(
      "Marihuana" = 0.8214285714,
      "Alcohol" = 0.6428571429,
      "Cocaína" = 0.2857142857,
      "Sedantes" = 0.1071428571,
      "Anfetaminas" = 0.0714285714,
      "Éxtasis/MDMA" = 0.0714285714,
      "Opioides" = 0.0357142857
    ),
    "Heroína" = list(
      "Cocaína" = 0.6666666667,
      "Marihuana" = 0.3888888889,
      "Alcohol" = 0.2777777778,
      "Anfetaminas" = 0.1666666667,
      "Sedantes" = 0.1666666667,
      "Crack" = 0.1111111111,
      "Éxtasis/MDMA" = 0.1111111111,
      "Alucinógenos" = 0.0555555556,
      "Lsd" = 0.0555555556,
      "Metanfetaminas y Otros Derivados" = 0.0555555556,
      "Opioides" = 0.0555555556,
      "Pasta Base" = 0.0555555556
    ),
    "Fenilciclidina" = list(
      "Alcohol" = 0.6000000000,
      "Cocaína" = 0.6000000000,
      "Heroína" = 0.4000000000,
      "Éxtasis/MDMA" = 0.2000000000
    ),
    "Metadona" = list(
      "Sedantes" = 0.2000000000
    ),
    "Esteroides Anabólicos" = list(
      "Alcohol" = 1.0000000000,
      "Cocaína" = 1.0000000000
    )
  )
  
  # Para sustancias sin probabilidades condicionales específicas
  dist_general_otras <- c(
    "Alcohol" = 0.35,
    "Marihuana" = 0.25,
    "Cocaína" = 0.20,
    "Pasta Base" = 0.10,
    "Sedantes" = 0.05,
    "Otros" = 0.02,
    "Opioides" = 0.015,
    "Estimulantes" = 0.01,
    "Hipnóticos" = 0.005
  )
  
  # ===============================================
  # GENERACIÓN DE DATOS
  # ===============================================
  
  # Inicializar dataframe
  data_sim <- data.frame(
    HASH_KEY = paste0("SIM_", sprintf("%08d", 1:n_pacientes)),
    sustancia_principal = NA_character_,
    otras_sustancias_no1 = NA_character_,
    otras_sustancias_no2 = NA_character_,
    otras_sustancias_no3 = NA_character_,
    stringsAsFactors = FALSE
  )
  
  # Asignar sustancias principales según distribución exacta
  data_sim$sustancia_principal <- sample(names(dist_principal),
                                         size = n_pacientes,
                                         prob = dist_principal,
                                         replace = TRUE)
  
  # Para cada paciente, determinar policonsumo
  for(i in 1:n_pacientes) {
    sust_prin <- data_sim$sustancia_principal[i]
    
    # Obtener probabilidad de policonsumo para esta sustancia
    prob_poli <- ifelse(sust_prin %in% names(prob_policonsumo),
                        prob_policonsumo[sust_prin],
                        0.70)
    
    # Determinar si tiene policonsumo
    if(runif(1) < prob_poli) {
      # Tiene policonsumo - determinar cuántas sustancias en total (2, 3 o 4)
      probs_n <- dist_n_sustancias[2:4]
      probs_n <- as.numeric(probs_n) / sum(as.numeric(probs_n))
      
      n_total <- sample(2:4, size = 1, prob = probs_n)
      n_adicionales <- n_total - 1
      
      # Obtener probabilidades condicionales para esta sustancia principal
      if(sust_prin %in% names(prob_condicionales)) {
        probs <- prob_condicionales[[sust_prin]]
      } else {
        probs <- dist_general_otras
      }
      
      # Verificar que hay sustancias disponibles
      if(length(probs) > 0) {
        # Convertir a vector con nombres
        probs_vec <- unlist(probs)
        
        # Eliminar la sustancia principal si aparece en las opciones
        probs_vec <- probs_vec[names(probs_vec) != sust_prin]
        
        # Eliminar sustancias con probabilidad 0 o NA
        probs_vec <- probs_vec[!is.na(probs_vec) & probs_vec > 0]
        
        if(length(probs_vec) >= n_adicionales) {
          # Seleccionar sustancias adicionales según probabilidades
          sust_seleccionadas <- sample(names(probs_vec),
                                      size = n_adicionales,
                                      prob = probs_vec,
                                      replace = FALSE)
          
          # Asignar a las columnas correspondientes
          if(length(sust_seleccionadas) >= 1) {
            data_sim$otras_sustancias_no1[i] <- sust_seleccionadas[1]
          }
          if(length(sust_seleccionadas) >= 2) {
            data_sim$otras_sustancias_no2[i] <- sust_seleccionadas[2]
          }
          if(length(sust_seleccionadas) >= 3) {
            data_sim$otras_sustancias_no3[i] <- sust_seleccionadas[3]
          }
        }
      }
    }
  }
  
  return(data_sim)
}

# ===============================================
# DETECCIÓN Y CARGA/SIMULACIÓN DE DATOS
# ===============================================

# Intentar cargar datos reales
data_path <- paste0(gsub("/docs", "", getwd()), "/data/CONS_C1_2010_22_CLEAN.rds")

if(file.exists(data_path)) {
  data <- readRDS(data_path)
} else {
  # Simular datos con propiedades idénticas a los reales
  n_sim <- 100000  # Se puede ajustar este número (original: 223055)
  data <- simular_datos_chile(n_pacientes = n_sim, seed = 2024)
}

# Función de simplificación
simplify_substance_names <- function(x) {
  x <- as.character(x)
  x[str_detect(x, "^Sin ")] <- NA
  x <- str_replace(x, "^Sedantes:.*", "Sedantes")
  x <- str_replace(x, "^Hipnóticos:.*", "Hipnóticos")
  x <- str_replace(x, "^Inhalables:.*", "Inhalables")
  x <- str_replace(x, "^Otros Opioides.*", "Opioides")
  x <- str_replace(x, "^Otros Estimulantes.*", "Estimulantes")
  x <- str_replace(x, "^Otros Alucinógenos.*", "Alucinógenos")
  x <- str_replace(x, "^Éxtasis.*", "Éxtasis/MDMA")
  return(x)
}

# Columnas de sustancias
cols_sustancias <- c("sustancia_principal", "otras_sustancias_no1", 
                     "otras_sustancias_no2", "otras_sustancias_no3")

# Aplicar limpieza
data_network <- data %>%
  select(HASH_KEY, all_of(cols_sustancias)) %>%
  mutate(across(all_of(cols_sustancias), simplify_substance_names)) %>%
  filter(!is.na(sustancia_principal))

1.2 INTRODUCCIÓN

Los trastornos por uso de sustancias (TUS) constituyen una de las principales causas de carga de enfermedad y mortalidad evitable a escala global. En 2016 se calculó que más de 100 millones de personas sufrían trastorno por consumo de alcohol y decenas de millones presentaban dependencia de opioides, cannabis o cocaína1. La frecuente comorbilidad psiquiátrica, depresión, trastornos de ansiedad, psicosis o trastornos de personalidad multiplica la severidad clínica y los costes sociosanitarios2. Estudios hospitalarios europeos y norteamericanos muestran que alrededor del 20 % de las admisiones psiquiátricas corresponden a pacientes de sexo femenino con diagnóstico dual, fenómeno que favorece re-ingresos y estancias prolongadas3.

En Chile, las encuestas nacionales sitúan la prevalencia de abuso o dependencia de sustancias entre el 11 % y el 20 %, una de las más elevadas de Latinoamérica. Los registros hospitalarios concuerdan con las cifras internacionales: alrededor de una quinta parte de los internados en psiquiatría presenta un TCS como diagnóstico primario o secundario. Esta convergencia evidencia que la hospitalización psiquiátrica es un desenlace clínico crítico en la trayectoria de las adicciones, razón por la cual identificar sus factores determinantes resulta esencial para planificar intervenciones preventivas, asignar recursos y reducir la carga asistencial2,46.

1.3 DESCRIPCIÓN GENERAL DE LOS DATOS

Este es un estudio de cohorte retrospectiva de pacientes adultos en tratamiento por consumo de sustancias, con datos otorgados por el Servicio Nacional para la Prevención y Rehabilitación del Consumo de Drogas y Alcohol de Chile (SENDA) en convenio con el núcleo milenio de ánalisis de políticas públicas de drogas (nDP). La cohorte se construyó vinculando los registros administrativos de los pacientes (n = 223,061 episodios de tratamiento entre 97,698 personas en las 16 regiones del país).

Estos datos incluyen múltiples variables relacionadas al consumo y tratamiento rehabilitador de drogas. Entre estas variables esta la sustancia principal por la cual se trató al paciente y sustancias secundarias (alcohol, pasta base de cocaína, cocaína, marihuana, depresores del SNC u otras sustancias. Tambien está presente el número de reingresos a tratamiento (retratamientos, categorizados en 0, 1, 2, 3 o más reingresos), el tipo de plan de tratamiento (ambulatorio vs. residencial) y el historial clínico de salud mental de los pacientes.

El registro de pacientes en tratamiento se realizó en una plataforma electrónica denominada SISTRAT, que contenía información sociodemográfica, datos sobre el estado de salud y patrones de consumo de sustancias, entre otras variables, además de información sobre el propio tratamiento (p. ej., fecha de ingreso, egreso, tipo de tratamiento). Las base de datos se vincularon de forma determinista mediante un hash de 64 caracteres resultante del cifrado (con un algoritmo SHA-256) del número de identificación único de cada persona.

Código
info_data <- data.frame(
  Característica = c("Total de registros",
                     "Registros con sustancia principal",
                     "Período de análisis",
                     "Variables de sustancias",
                     "Número de sustancias únicas",
                     "Promedio de sustancias por paciente",
                     "Mediana de sustancias por paciente",
                     "Desviación estándar"),
  Valor = c(format(nrow(data), big.mark = ","),
            format(nrow(data_network), big.mark = ","),
            "2010-2022",
            as.character(length(cols_sustancias)),
            as.character(n_distinct(unlist(data_network[cols_sustancias]), na.rm = TRUE)),
            round(mean(rowSums(!is.na(data_network[cols_sustancias]))), 2),
            median(rowSums(!is.na(data_network[cols_sustancias]))),
            round(sd(rowSums(!is.na(data_network[cols_sustancias]))), 2))
)

info_data %>%
  kable(format = "html", 
        col.names = c("Característica", "Valor"),
        align = c("l", "r"),
        row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE)
Características generales de la base de datos
Característica Valor
Total de registros 223,061
Registros con sustancia principal 223,055
Período de análisis 2010-2022
Variables de sustancias 4
Número de sustancias únicas 22
Promedio de sustancias por paciente 2.2
Mediana de sustancias por paciente 2
Desviación estándar 0.99

2 RED DE CO-OCURRENCIA

2.1 Construcción de la Red

La red de co-ocurrencia es pertinente para identificar patrones de policonsumo y asociaciones entre sustancias que tienden a ser usadas conjuntamente. Esta representación revela combinaciones frecuentes y comunidades de sustancias; en particular, ayuda a localizar “nodos puente” que conectan módulos y que suelen ser dianas clínicas de alto impacto para el cribado y la intervención integrada79. Además, varios emparejamientos visibles en la red coinciden con riesgos clínicos conocidos: el uso concurrente de opioides y benzodiacepinas se asocia con un aumento de 2–5 veces del riesgo de sobredosis (mayor en los primeros 90 días de co‑exposición), y más del 90% de las muertes con benzodiacepinas co‑involucran opioides1012; la co‑ingesta de alcohol y cocaína genera cocaetileno, metabolito más cardiotóxico, por lo que constituye una combinación especialmente peligrosa13; y la co‑involucración de estimulantes con opioides sintéticos (p. ej., fentanilo) ha aumentado de forma marcada en los últimos años, subrayando la necesidad de estrategias específicas para policonsumo14. Metodológicamente, la transformación logarítmica de pesos resulta adecuada en redes con distribuciones de cola pesada, permitiendo visualizar de forma equilibrada tanto las asociaciones muy frecuentes como las moderadas15. Por su parte, filtrar el 60% superior de co‑ocurrencias es coherente con los enfoques de “backbone” en redes ponderadas, que eliminan enlaces débiles para resaltar la estructura multiescala estadísticamente significativa y focalizar los patrones robustos16.

Código
# Crear matriz de co-ocurrencia
create_cooccurrence_matrix <- function(df) {
  all_substances <- unique(unlist(df[cols_sustancias]))
  all_substances <- all_substances[!is.na(all_substances)]
  all_substances <- sort(all_substances)
  
  n_sust <- length(all_substances)
  co_matrix <- matrix(0, nrow = n_sust, ncol = n_sust,
                     dimnames = list(all_substances, all_substances))
  
  for (i in 1:nrow(df)) {
    row_substances <- unlist(df[i, cols_sustancias])
    row_substances <- row_substances[!is.na(row_substances)]
    
    if (length(row_substances) > 1) {
      for (j in 1:(length(row_substances)-1)) {
        for (k in (j+1):length(row_substances)) {
          sust1 <- row_substances[j]
          sust2 <- row_substances[k]
          
          co_matrix[sust1, sust2] <- co_matrix[sust1, sust2] + 1
          co_matrix[sust2, sust1] <- co_matrix[sust2, sust1] + 1
        }
      }
    }
  }
  
  return(co_matrix)
}

co_matrix <- create_cooccurrence_matrix(data_network)

# IMPORTANTE: Solo el 60% de las co-ocurrencias más frecuentes
# Obtener todos los valores únicos de co-ocurrencia (excluyendo diagonal)
co_values <- co_matrix[upper.tri(co_matrix)]
co_values <- co_values[co_values > 0]

# Calcular el percentil 40 (para mantener el 60% superior)
threshold_value <- quantile(co_values, probs = 0.40)

# Aplicar umbral a la matriz
co_matrix_filtered <- co_matrix
co_matrix_filtered[co_matrix < threshold_value] <- 0

# Crear grafo con pesos logarítmicos (ESCALA LOGARÍTMICA)
g_full_co <- graph_from_adjacency_matrix(
  co_matrix_filtered,
  mode = "undirected",
  weighted = TRUE,
  diag = FALSE
)

# Aplicar transformación logarítmica a los pesos
E(g_full_co)$weight_original <- E(g_full_co)$weight
E(g_full_co)$weight <- log1p(E(g_full_co)$weight)  # ESCALA LOGARÍTMICA: log(1 + weight)

# Eliminar vértices aislados
g_co <- delete_vertices(g_full_co, degree(g_full_co) == 0)

# Detectar comunidades
communities_co <- cluster_louvain(g_co)
V(g_co)$community <- membership(communities_co)

# Información sobre el filtrado
n_edges_original <- sum(co_matrix > 0) / 2
n_edges_filtered <- sum(co_matrix_filtered > 0) / 2
percent_retained <- round(n_edges_filtered / n_edges_original * 100, 1)

2.2 Métricas de la Red de Co-ocurrencia

Código
# Calcular métricas de centralidad
centrality_metrics_co <- data.frame(
  Sustancia = V(g_co)$name,
  Grado = degree(g_co),
  Fuerza = round(strength(g_co), 2),
  Intermediación = round(betweenness(g_co, normalized = TRUE), 3),
  Cercanía = round(closeness(g_co, normalized = TRUE), 3),
  Eigenvector = round(eigen_centrality(g_co)$vector, 3),
  PageRank = round(page_rank(g_co)$vector, 4),
  Comunidad = V(g_co)$community
) %>%
  arrange(desc(Fuerza))

# Propiedades estructurales
network_props_co <- data.frame(
  Propiedad = c("Número de nodos", "Número de enlaces", 
                "Enlaces filtrados (60% superior)", "Densidad",
                "Diámetro", "Distancia media", "Clustering global",
                "Modularidad", "Asortatividad"),
  Valor = c(vcount(g_co), ecount(g_co), 
            paste0(n_edges_filtered, " de ", n_edges_original, " (", percent_retained, "%)"),
            round(edge_density(g_co), 4),
            diameter(g_co, weights = NA),
            round(mean_distance(g_co, weights = NA), 2),
            round(transitivity(g_co, type = "global"), 3),
            round(modularity(communities_co), 3),
            round(assortativity_degree(g_co), 3))
)

# Mostrar top 10 sustancias
centrality_metrics_co %>%
  head(10) %>%
  select(Sustancia, Grado, Fuerza, Intermediación, Cercanía, Eigenvector, PageRank) %>%
  kable(format = "html", row.names = F) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Sustancia Grado Fuerza Intermediación Cercanía Eigenvector PageRank
Alcohol 20 131.26 0.021 0.174 1.000 0.1106
Cocaína 20 125.12 0.137 0.180 0.960 0.1056
Marihuana 19 124.06 0.042 0.170 0.971 0.1024
Pasta Base 18 112.01 0.116 0.165 0.918 0.0917
Sedantes 15 79.09 0.000 0.163 0.720 0.0651
Opioides 15 60.63 0.116 0.193 0.557 0.0518
Otros 12 59.05 0.011 0.159 0.590 0.0499
Éxtasis/MDMA 14 56.96 0.026 0.191 0.526 0.0488
Anfetaminas 12 52.92 0.016 0.171 0.535 0.0453
Lsd 11 45.77 0.011 0.184 0.455 0.0403
Código
# Mostrar propiedades
network_props_co %>%
  kable(format = "html", align = c("l", "r")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Propiedad Valor
Número de nodos 21
Número de enlaces 114
Enlaces filtrados (60% superior) 116.5 de 185 (63%)
Densidad 0.5429
Diámetro 2
Distancia media 1.46
Clustering global 0.699
Modularidad 0.018
Asortatividad -0.488

Métricas principales de la red de co-ocurrencia (60% de co-ocurrencias más frecuentes)

2.3 Visualización Estática

Código
g_co_tidy <- as_tbl_graph(g_co)

# Calcular tamaños proporcionales a la fuerza (con pesos logarítmicos)
node_strength <- strength(g_co)
node_size <- sqrt(node_strength) / max(sqrt(node_strength)) * 20 + 5

set.seed(42)
ggraph(g_co_tidy, layout = 'fr', weights = weight) + 
  geom_edge_link(aes(width = weight, alpha = weight),
                 color = "gray40",
                 show.legend = FALSE) +
  geom_node_point(aes(color = factor(V(g_co)$community)),
                  size = node_size,
                  alpha = 0.9) +
  geom_node_text(aes(label = name),
                 size = node_size/4,
                 repel = TRUE,
                 force = 3,
                 segment.size = 0.2,
                 segment.alpha = 0.5,
                 max.overlaps = Inf, 
                 show.legend = FALSE) +
  scale_edge_width_continuous(range = c(0.2, 3)) +
  scale_edge_alpha_continuous(range = c(0.2, 0.7)) +
  scale_color_brewer(palette = "Set1",
                     name = "Comunidad") +
  labs(title = "Red de Co-ocurrencia de Sustancias",
       subtitle = paste0("Escala logarítmica | 60% de co-ocurrencias más frecuentes (", 
                        n_edges_filtered, " de ", n_edges_original, " enlaces)")) +
  theme_void() +
  theme(plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        plot.subtitle = element_text(size = 12, hjust = 0.5),
        legend.position = "right")

Red de co-ocurrencia de sustancias (escala logarítmica, 60% superior de co-ocurrencias)

2.4 Red Interactiva Co-ocurrencia

Código
# Calcular layout con mayor separación
set.seed(123)
coords_co <- layout_with_fr(g_co, weights = E(g_co)$weight)
coords_co <- coords_co * 300

# Calcular tamaños proporcionales
node_strength <- strength(g_co)
node_sizes <- (node_strength / max(node_strength)) * 50 + 10

# Calcular k-core
coreness_co <- coreness(g_co)

# Preparar datos para visNetwork
nodes_co <- data.frame(
  id = V(g_co)$name,
  label = V(g_co)$name,
  value = node_sizes,
  group = V(g_co)$community,
  x = coords_co[,1],
  y = coords_co[,2],
  title = paste0(
    "<b>", V(g_co)$name, "</b><br>",
    "Fuerza (log): ", format(round(strength(g_co), 2), big.mark = ","), "<br>",
    "Co-ocurrencias originales: ", format(round(sum(exp(E(g_co)$weight[E(g_co)$weight > 0]) - 1)), big.mark = ","), "<br>",
    "Conexiones: ", degree(g_co), "<br>",
    "Comunidad: ", V(g_co)$community, "<br>",
    "K-Core: ", coreness_co
  ),
  font.size = node_sizes * 0.8,  # Tamaño de letra proporcional al nodo
  font.face = "Arial",  # Fuente Arial
  physics = TRUE
)

# Ajustar posiciones de nodos principales
if ("Alcohol" %in% nodes_co$id) {
  idx <- which(nodes_co$id == "Alcohol")
  nodes_co[idx, c("x", "y")] <- c(-400, 0)
  nodes_co[idx, "fixed"] <- TRUE
}
if ("Marihuana" %in% nodes_co$id) {
  idx <- which(nodes_co$id == "Marihuana")
  nodes_co[idx, c("x", "y")] <- c(400, 0)
  nodes_co[idx, "fixed"] <- TRUE
}
if ("Pasta Base" %in% nodes_co$id) {
  idx <- which(nodes_co$id == "Pasta Base")
  nodes_co[idx, c("x", "y")] <- c(0, 400)
  nodes_co[idx, "fixed"] <- TRUE
}
if ("Cocaína" %in% nodes_co$id) {
  idx <- which(nodes_co$id == "Cocaína")
  nodes_co[idx, c("x", "y")] <- c(0, -400)
  nodes_co[idx, "fixed"] <- TRUE
}

# Preparar edges con grosor proporcional (usando pesos originales para el tooltip)
edges_df_co <- as_data_frame(g_co, what = "edges")
edges_co <- data.frame(
  from = edges_df_co$from,
  to = edges_df_co$to,
  value = edges_df_co$weight / 5,  # Escalar peso logarítmico para visualización
  title = paste0(
    edges_df_co$from, " ↔ ", edges_df_co$to, "<br>",
    "Co-ocurrencias: ", format(round(exp(edges_df_co$weight) - 1), big.mark = ","), "<br>",
    "Peso logarítmico: ", round(edges_df_co$weight, 2)
  ),
  color = "rgba(150,150,150,0.3)",
  highlight = "rgba(255,100,100,0.8)"
)

# Crear visualización interactiva
visNetwork(nodes_co, edges_co, height = "700px", width = "100%") %>%
  visOptions(
    highlightNearest = list(
      enabled = TRUE,
      hover = TRUE,
      degree = 1,
      labelOnly = FALSE
    ),
    nodesIdSelection = list(
      enabled = TRUE,
      style = 'width: 200px; height: 30px;'
    ),
    selectedBy = list(
      variable = "group",
      style = 'width: 200px; height: 30px;',
      main = "Seleccionar Comunidad"
    )
  ) %>%
  visPhysics(
    enabled = TRUE,
    stabilization = list(
      enabled = TRUE,
      iterations = 500,
      updateInterval = 50
    ),
    barnesHut = list(
      gravitationalConstant = -8000,
      centralGravity = 0.1,
      springLength = 300,
      springConstant = 0.001,
      damping = 0.5,
      avoidOverlap = 1
    ),
    maxVelocity = 50,
    minVelocity = 0.75,
    solver = "barnesHut"
  ) %>%
  visLayout(
    randomSeed = 123,
    improvedLayout = TRUE
  ) %>%
  visInteraction(
    navigationButtons = TRUE,
    dragNodes = TRUE,
    dragView = TRUE,
    zoomView = TRUE,
    hover = TRUE,
    tooltipDelay = 0,
    zoomSpeed = 0.5
  ) %>%
  visEdges(
    smooth = list(
      enabled = TRUE,
      type = "continuous",
      roundness = 0.5
    ),
    width = 2,
    physics = TRUE,
    scaling = list(
      min = 0.5,
      max = 5
    )
  ) %>%
  visNodes(
    shape = "dot",
    scaling = list(
      min = 10,
      max = 60,
      label = list(
        enabled = TRUE,
        min = 10,
        max = 40,
        drawThreshold = 0,
        maxVisible = 50
      )
    ),
    font = list(
      face = "Arial, sans-serif",
      align = "center"
    ),
    borderWidth = 2,
    borderWidthSelected = 4
  ) %>%
  visLegend(
    enabled = TRUE,
    position = "right",
    width = 0.15,
    useGroups = TRUE,
    main = "Comunidades"
  ) %>%
  visConfigure(enabled = FALSE)

Red interactiva de co-ocurrencia (zoom y arrastre habilitados)

3 RED BIPARTITA PACIENTE-SUSTANCIA

3.1 Construcción de la Red

La red bipartita es útil y que representa simultáneamente a los pacientes y a las sustancias, preservando la heterogeneidad individual y, a la vez, permitiendo ver la prevalencia relativa de cada sustancia a través del grado de sus nodos17,18. Al mantener la estructura de dos modos, se evita la pérdida de información y las distorsiones que introducen las proyecciones unipartitas y se pueden aplicar métricas específicas de bipartitos para capturar patrones de alto riesgo19. Además, medidas como la betweenness centrality identifican sustancias que actúan como “puentes” entre grupos de consumidores, clave para detectar rutas de transición en el policonsumo20. Para descubrir perfiles y comunidades de policonsumo se dispone de modularidad bipartita y variantes para redes ponderadas, validadas en la literatura metodológica21,22. La utilidad clínica de los bipartitos está documentada en contextos de salud: redes paciente‑prescriptor y paciente‑hospital han permitido detectar conductas de búsqueda de fármacos y diferenciar patrones según sustancia, evidenciando nodos (sustancias o instituciones) estratégicos para intervención19,Kim2024?. Finalmente, para equilibrar legibilidad y representatividad de la visualización, el muestreo aleatorio simple de 2.500 pacientes facilita la visulización sin introducir sesgos y, usado con fines gráficos, preserva propiedades globales (distribuciones y relaciones clave) de forma aceptable para el análisis exploratorio23,24.

Código
# MUESTREO ALEATORIO SIMPLE DE 2,500 SUJETOS DEL TOTAL
# No estratificado por sustancia - cada paciente tiene igual probabilidad de ser seleccionado
set.seed(42)
n_sample <- 10000
sample_patients <- sample(unique(data_network$HASH_KEY), 
                         min(n_sample, length(unique(data_network$HASH_KEY))))
data_network_sampled <- data_network %>%
  filter(HASH_KEY %in% sample_patients)

# Crear edge list paciente-sustancia con datos muestreados
create_bipartite_edgelist <- function(data) {
  edgelist <- data.frame()
  
  for (col in cols_sustancias) {
    temp <- data %>%
      select(HASH_KEY, sustancia = !!sym(col)) %>%
      filter(!is.na(sustancia)) %>%
      mutate(weight = ifelse(col == "sustancia_principal", 2, 1))
    
    edgelist <- rbind(edgelist, temp)
  }
  
  edgelist <- edgelist %>%
    group_by(HASH_KEY, sustancia) %>%
    summarise(weight = sum(weight), .groups = 'drop')
  
  return(edgelist)
}

bipartite_edges <- create_bipartite_edgelist(data_network_sampled)

# Crear grafo bipartito
g_bipartite <- graph_from_data_frame(bipartite_edges, directed = FALSE)
V(g_bipartite)$type <- V(g_bipartite)$name %in% bipartite_edges$sustancia
V(g_bipartite)$node_type <- ifelse(V(g_bipartite)$type, "Sustancia", "Paciente")

# Detectar comunidades
communities_bi <- cluster_louvain(g_bipartite)
V(g_bipartite)$community <- membership(communities_bi)

# Estadísticas básicas
n_pacientes <- sum(!V(g_bipartite)$type)
n_sustancias <- sum(V(g_bipartite)$type)
n_conexiones <- ecount(g_bipartite)

3.2 Métricas de la Red Bipartita

Código
# Calcular grado por tipo
degree_all_bi <- degree(g_bipartite)
degree_pacientes <- degree_all_bi[!V(g_bipartite)$type]
degree_sustancias <- degree_all_bi[V(g_bipartite)$type]

# Propiedades estructurales
network_props_bi <- data.frame(
  Propiedad = c("Número de pacientes (muestra)", "Número de sustancias",
                "Total de conexiones", "Densidad",
                "Conexiones promedio por paciente",
                "Conexiones promedio por sustancia",
                "Componentes conectados", "Modularidad"),
  Valor = c(format(n_pacientes, big.mark = ","),
            n_sustancias,
            format(n_conexiones, big.mark = ","),
            round(n_conexiones / (n_pacientes * n_sustancias), 6),
            round(mean(degree_pacientes), 2),
            round(mean(degree_sustancias), 0),
            components(g_bipartite)$no,
            round(modularity(communities_bi), 3))
)

network_props_bi %>%
  kable(format = "html", align = c("l", "r")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Métricas principales de la red bipartita (muestra aleatoria simple de 10000 pacientes)
Propiedad Valor
Número de pacientes (muestra) 10,000
Número de sustancias 21
Total de conexiones 23,735
Densidad 0.113024
Conexiones promedio por paciente 2.37
Conexiones promedio por sustancia 1130
Componentes conectados 1
Modularidad 0.341

3.3 Visualización Estática Bipartita

Código
# Cargar ggrepel
require(ggrepel)

# Usar el grafo muestreado
g_bi_viz <- g_bipartite

# Calcular grados para TODOS los nodos
degree_bi_viz <- degree(g_bi_viz)

# Calcular tamaños proporcionales al grado
node_sizes_bi <- numeric(vcount(g_bi_viz))

# Para sustancias: escala más grande
substance_nodes <- which(V(g_bi_viz)$type)
if(length(substance_nodes) > 0) {
  substance_degrees <- degree_bi_viz[substance_nodes]
  node_sizes_bi[substance_nodes] <- 15 + (substance_degrees / max(substance_degrees)) * 30
}

# Para pacientes: escala más pequeña
patient_nodes <- which(!V(g_bi_viz)$type)
if(length(patient_nodes) > 0) {
  patient_degrees <- degree_bi_viz[patient_nodes]
  node_sizes_bi[patient_nodes] <- 0.5 + (patient_degrees / max(patient_degrees)) * 5
}

# ===== OPCIÓN 1: LAYOUT CIRCULAR/RADIAL =====
create_radial_layout <- function() {
  layout_matrix <- matrix(0, vcount(g_bi_viz), 2)
  
  # Pacientes en círculo interior (radio menor)
  if(length(patient_nodes) > 0) {
    angles_patients <- seq(0, 2*pi, length.out = length(patient_nodes) + 1)[1:length(patient_nodes)]
    radius_patients <- 3  # Radio interior
    
    # Variar ligeramente el radio según el grado para crear textura
    radius_variation <- 0.3 * (degree_bi_viz[patient_nodes] / max(degree_bi_viz[patient_nodes]))
    
    layout_matrix[patient_nodes, 1] <- (radius_patients + radius_variation) * cos(angles_patients)
    layout_matrix[patient_nodes, 2] <- (radius_patients + radius_variation) * sin(angles_patients)
  }
  
  # Sustancias en círculo exterior (radio mayor)
  if(length(substance_nodes) > 0) {
    # Ordenar por grado para mejor distribución
    substance_order <- substance_nodes[order(degree_bi_viz[substance_nodes], decreasing = TRUE)]
    
    # Distribuir ángulos con peso según el grado
    weights <- degree_bi_viz[substance_order] + 10
    cumulative_weights <- cumsum(weights)
    angles_substances <- (cumulative_weights / sum(weights)) * 2 * pi
    
    radius_substances <- 5  # Radio exterior
    
    for(i in seq_along(substance_order)) {
      layout_matrix[substance_order[i], 1] <- radius_substances * cos(angles_substances[i])
      layout_matrix[substance_order[i], 2] <- radius_substances * sin(angles_substances[i])
    }
  }
  
  return(layout_matrix)
}

# ===== OPCIÓN 2: LAYOUT ARCO =====
create_arc_layout <- function() {
  layout_matrix <- matrix(0, vcount(g_bi_viz), 2)
  
  # Pacientes en arco inferior
  if(length(patient_nodes) > 0) {
    # Crear un arco de -120 a 120 grados
    angles <- seq(-2*pi/3, 2*pi/3, length.out = length(patient_nodes))
    radius <- 4
    
    # Agregar variación en radio según grado
    radius_var <- radius + 0.5 * (degree_bi_viz[patient_nodes] / max(degree_bi_viz[patient_nodes]))
    
    layout_matrix[patient_nodes, 1] <- radius_var * cos(angles)
    layout_matrix[patient_nodes, 2] <- radius_var * sin(angles) - 2  # Desplazar hacia abajo
  }
  
  # Sustancias en arco superior
  if(length(substance_nodes) > 0) {
    substance_order <- substance_nodes[order(degree_bi_viz[substance_nodes], decreasing = TRUE)]
    
    # Distribuir con peso en un arco superior más pequeño
    weights <- sqrt(degree_bi_viz[substance_order]) + 1
    positions <- cumsum(c(0, weights[-length(weights)]))
    positions <- (positions / max(positions)) * pi - pi/2  # Arco de -90 a 90 grados
    
    radius <- 6
    
    for(i in seq_along(substance_order)) {
      layout_matrix[substance_order[i], 1] <- radius * cos(positions[i])
      layout_matrix[substance_order[i], 2] <- radius * sin(positions[i]) + 3  # Desplazar hacia arriba
    }
  }
  
  return(layout_matrix)
}

# ===== OPCIÓN 3: LAYOUT GRID EXPANDIDO =====
create_expanded_grid_layout <- function() {
  layout_matrix <- matrix(0, vcount(g_bi_viz), 2)
  
  # Pacientes en grid expandido con jitter
  if(length(patient_nodes) > 0) {
    n_patients <- length(patient_nodes)
    # Crear una cuadrícula más dispersa
    n_cols <- ceiling(sqrt(n_patients * 2))  # Más columnas para mayor dispersión
    n_rows <- ceiling(n_patients / n_cols)
    
    # Ordenar por grado para agrupar similares
    patient_order <- patient_nodes[order(degree_bi_viz[patient_nodes], decreasing = TRUE)]
    
    idx <- 1
    for(i in 1:n_rows) {
      for(j in 1:n_cols) {
        if(idx <= n_patients) {
          # Agregar jitter proporcional al grado
          jitter_strength <- 0.2 * (degree_bi_viz[patient_order[idx]] / max(degree_bi_viz[patient_nodes]))
          layout_matrix[patient_order[idx], 1] <- (j - n_cols/2) * 0.3 + runif(1, -jitter_strength, jitter_strength)
          layout_matrix[patient_order[idx], 2] <- (i - n_rows/2) * 0.3 - 2  # Abajo
          idx <- idx + 1
        }
      }
    }
  }
  
  # Sustancias en arco superior expandido
  if(length(substance_nodes) > 0) {
    substance_order <- substance_nodes[order(degree_bi_viz[substance_nodes], decreasing = TRUE)]
    
    # Distribuir en dos filas escalonadas
    n_sust <- length(substance_order)
    half <- ceiling(n_sust / 2)
    
    # Primera fila (sustancias más importantes)
    for(i in 1:min(half, n_sust)) {
      x_pos <- (i - half/2 - 0.5) * 1.2
      layout_matrix[substance_order[i], 1] <- x_pos
      layout_matrix[substance_order[i], 2] <- 3
    }
    
    # Segunda fila (si hay suficientes sustancias)
    if(n_sust > half) {
      for(i in (half+1):n_sust) {
        x_pos <- ((i-half) - (n_sust-half)/2 - 0.5) * 1.2
        layout_matrix[substance_order[i], 1] <- x_pos
        layout_matrix[substance_order[i], 2] <- 2.2
      }
    }
  }
  
  return(layout_matrix)
}

# SELECCIONAR EL LAYOUT (cambiar según preferencia)
# Opción 1: Radial/Circular
# layout_matrix <- create_radial_layout()

# Opción 2: Arco
# layout_matrix <- create_arc_layout()

# Opción 3: Grid Expandido (por defecto)
set.seed(42)  # Para reproducibilidad del jitter
layout_matrix <- create_expanded_grid_layout()

# Crear el grafo tidy
g_bi_tidy <- as_tbl_graph(g_bi_viz) |>
  activate(nodes) |>
  mutate(
    node_type = ifelse(type, "Sustancia", "Paciente"),
    node_size_plot = node_sizes_bi,
    node_degree = degree_bi_viz,
    label = ifelse(type, name, NA_character_),
    node_color = ifelse(type, "#e74c3c", "#3498db")
  )

# Crear la visualización
p <- ggraph(g_bi_tidy, layout = layout_matrix) +
  # Edges con curvatura para mejor visualización
  geom_edge_link(aes(alpha = weight),
                 color = "gray75", 
                 width = 0.1,
                 show.legend = FALSE) +
  
  # Nodos
  geom_node_point(aes(size = node_size_plot, 
                      color = node_type,
                      alpha = node_type),
                  show.legend = TRUE) +
  
  # Etiquetas con ggrepel
  ggrepel::geom_text_repel(data = function(x) filter(x, !is.na(label)),
                            aes(x = x, y = y, label = label),
                            size = 3.5,
                            box.padding = 0.5,
                            point.padding = 0.3,
                            segment.color = "gray50",
                            segment.size = 0.3,
                            segment.alpha = 0.5,
                            max.overlaps = 50,
                            fontface = "bold",
                            force = 2,  # Mayor fuerza de repulsión
                            force_pull = 0.5,
                            max.time = 2,
                            max.iter = 10000) +
  
  scale_size_identity() +
  scale_color_manual(values = c("Paciente" = "#3498db", 
                                "Sustancia" = "#e74c3c"),
                     name = "Tipo de Nodo") +
  scale_alpha_manual(values = c("Paciente" = 0.5, 
                                "Sustancia" = 0.9),
                     guide = "none") +
  scale_edge_alpha_continuous(range = c(0.02, 0.15)) +
  
  # Ajustar límites para dar más espacio
  coord_cartesian(xlim = c(-8, 8), ylim = c(-5, 6)) +
  
  labs(title = "Red Bipartita Paciente-Sustancia",
       subtitle = paste0("Muestra: ",
                         format(n_pacientes, big.mark = ","),
                         " pacientes y ", n_sustancias, " sustancias | ",
                         format(n_conexiones, big.mark = ","), " conexiones"),
       caption = "El tamaño de cada nodo es proporcional a su número de conexiones") +
  theme_void() +
  theme(plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
        plot.subtitle = element_text(size = 13, hjust = 0.5, color = "gray40"),
        plot.caption = element_text(size = 10, hjust = 0.5, color = "gray50", 
                                   margin = margin(t = 10)),
        plot.background = element_rect(fill = "white", color = NA),
        legend.position = "bottom",
        legend.title = element_text(size = 11),
        legend.text = element_text(size = 10))

print(p)

Red bipartita paciente-sustancia (muestra aleatoria simple de 2,500 pacientes)

3.4 Red Interactiva Bipartita

Código
# Usar el grafo muestreado
g_bi_interactive <- g_bipartite

# Calcular grados
degree_interactive <- degree(g_bi_interactive)

# Preparar datos para visNetwork
nodes_bi_viz <- data.frame(
  id = V(g_bi_interactive)$name,
  label = ifelse(V(g_bi_interactive)$type, 
                 V(g_bi_interactive)$name, 
                 ""),
  group = V(g_bi_interactive)$node_type,
  value = ifelse(V(g_bi_interactive)$type,
                 degree_interactive * 2,
                 1),
  title = paste0(
    ifelse(V(g_bi_interactive)$type,
           paste0("<b>", V(g_bi_interactive)$name, "</b><br>",
                  "Total pacientes: ", format(degree_interactive, big.mark = ",")),
           paste0("Paciente ID: ", substr(V(g_bi_interactive)$name, 1, 8), "<br>",
                  "Sustancias: ", degree_interactive))
  ),
  shape = ifelse(V(g_bi_interactive)$type, "dot", "dot"),
  font.size = ifelse(V(g_bi_interactive)$type, 
                     sqrt(degree_interactive) * 3,  # Tamaño proporcional para sustancias
                     0),  # Sin etiquetas para pacientes
  font.face = "Arial",  # Fuente Arial
  physics = TRUE
)

# Layout bipartito manual mejorado
layout_bi <- matrix(0, nrow(nodes_bi_viz), 2)
substance_idx <- which(V(g_bi_interactive)$type)
patient_idx <- which(!V(g_bi_interactive)$type)

# SUSTANCIAS ARRIBA, con mayor separación
if(length(substance_idx) > 0) {
  sust_degrees <- degree_interactive[substance_idx]
  sust_order <- substance_idx[order(sust_degrees, decreasing = TRUE)]
  
  weights <- sqrt(sust_degrees[order(sust_degrees, decreasing = TRUE)])
  positions <- cumsum(c(0, weights[-length(weights)]))
  positions <- (positions / max(positions) - 0.5) * 1700  # mayor separación
  
  for(i in seq_along(sust_order)) {
    layout_bi[sust_order[i], 1] <- positions[i]
    layout_bi[sust_order[i], 2] <- -300  # Sustancias abajo
  }
}

# PACIENTES ABAJO, distribuidos con mayor separación
if(length(patient_idx) > 0) {
  n_patients <- length(patient_idx)
  
  if(n_patients > 1000) {
    n_rows <- ceiling(n_patients / 150)  # Reducido de 200 a 150 columnas para mayor dispersión
    for(i in seq_along(patient_idx)) {
      row <- ((i - 1) %/% 150) + 1
      col <- ((i - 1) %% 150) + 1
      layout_bi[patient_idx[i], 1] <- (col - 75) * 12  # Aumentado espaciado de 8 a 12
      layout_bi[patient_idx[i], 2] <- 300 + row * 60   # Aumentado espaciado vertical a 60
    }
  } else {
    layout_bi[patient_idx, 1] <- seq(-1200, 1200, length.out = n_patients)  # Aumentado rango
    layout_bi[patient_idx, 2] <- 300  # Pacientes arriba
  }
}

nodes_bi_viz$x <- layout_bi[,1]
nodes_bi_viz$y <- layout_bi[,2]

# Edges para visNetwork
edges_bi_viz <- as_data_frame(g_bi_interactive, what = "edges") %>%
  mutate(
    from = from,
    to = to,
    color = "rgba(150,150,150,0.15)",
    highlight = "rgba(255,107,107,0.8)",
    width = 0.5,
    smooth = FALSE
  )

# Crear red interactiva
visNetwork(nodes_bi_viz, edges_bi_viz, height = "800px", width = "100%") %>%
  visOptions(
    highlightNearest = list(
      enabled = TRUE,
      hover = TRUE,
      degree = 1,
      labelOnly = FALSE
    ),
    nodesIdSelection = list(
      enabled = TRUE,
      style = 'width: 200px; height: 30px;',
      main = "Buscar nodo"
    ),
    selectedBy = list(
      variable = "group",
      style = 'width: 200px; height: 30px;',
      main = "Filtrar por tipo"
    )
  ) %>%
  visPhysics(
    enabled = FALSE,
    stabilization = FALSE
  ) %>%
  visInteraction(
    navigationButtons = TRUE,
    dragNodes = TRUE,
    dragView = TRUE,
    zoomView = TRUE,
    hover = TRUE,
    tooltipDelay = 0,
    hideEdgesOnDrag = TRUE,
    hideNodesOnDrag = FALSE
  ) %>%
  visEdges(
    smooth = FALSE,
    physics = FALSE,
    hidden = FALSE,
    selectionWidth = 2
  ) %>%
  visNodes(
    scaling = list(
      min = 2,
      max = 60,
      label = list(
        enabled = TRUE,
        min = 8,
        max = 35,
        drawThreshold = 0,
        maxVisible = 30
      )
    ),
    font = list(
      face = "Arial, sans-serif",
      align = "center"
    ),
    borderWidth = 1,
    borderWidthSelected = 3
  ) %>%
  visGroups(
    groupname = "Sustancia",
    color = list(
      background = "#e74c3c",
      border = "#c0392b",
      highlight = list(
        background = "#ff6b6b",
        border = "#e74c3c"
      )
    ),
    shape = "dot",
    font = list(
      size = 16, 
      color = "black", 
      face = "Arial, sans-serif"
    )
  ) %>%
  visGroups(
    groupname = "Paciente",
    color = list(
      background = "#3498db",
      border = "#2980b9",
      highlight = list(
        background = "#5dade2",
        border = "#3498db"
      )
    ),
    shape = "dot",
    size = 3
  ) %>%
  visLegend(
    enabled = TRUE,
    position = "right",
    width = 0.1,
    useGroups = TRUE,
    main = "Tipo de Nodo"
  ) %>%
  visConfigure(enabled = FALSE)

Red bipartita interactiva (muestra aleatoria simple de 2,500 pacientes - zoom y arrastre habilitados)

4 RED BIPARTITA PROYECTADA

4.1 Construcción de la Proyección

Código
# NOTA: Para la proyección usamos TODOS los datos, no la muestra
# Crear edge list completo para proyección
bipartite_edges_full <- create_bipartite_edgelist(data_network)

# Crear proyección en sustancias
create_substance_projection <- function(edges) {
  substances <- unique(edges$sustancia)
  n_sust <- length(substances)
  
  adj_matrix <- matrix(0, n_sust, n_sust,
                      dimnames = list(substances, substances))
  
  for (patient in unique(edges$HASH_KEY)) {
    patient_subs <- edges$sustancia[edges$HASH_KEY == patient]
    if (length(patient_subs) > 1) {
      for (i in 1:(length(patient_subs)-1)) {
        for (j in (i+1):length(patient_subs)) {
          adj_matrix[patient_subs[i], patient_subs[j]] <- 
            adj_matrix[patient_subs[i], patient_subs[j]] + 1
          adj_matrix[patient_subs[j], patient_subs[i]] <- 
            adj_matrix[patient_subs[j], patient_subs[i]] + 1
        }
      }
    }
  }
  
  return(adj_matrix)
}

# Usar todos los datos para la proyección
adj_substances <- create_substance_projection(bipartite_edges_full)
g_substances <- graph_from_adjacency_matrix(adj_substances, 
                                           mode = "undirected", 
                                           weighted = TRUE)

# Detectar comunidades
communities_proj <- cluster_louvain(g_substances)
V(g_substances)$community <- membership(communities_proj)

4.2 Visualización Estática de Proyección

Código
g_proj_tidy <- as_tbl_graph(g_substances)

# Calcular tamaños proporcionales
node_strength_proj <- strength(g_substances)
node_size_proj <- sqrt(node_strength_proj) / max(sqrt(node_strength_proj)) * 20 + 5

set.seed(42)
ggraph(g_proj_tidy, layout = 'fr', weights = weight) + 
  geom_edge_link(aes(width = weight, alpha = weight),
                 color = "gray40",
                 show.legend = FALSE) +
  geom_node_point(aes(color = factor(V(g_substances)$community)),
                  size = node_size_proj,
                  alpha = 0.9) +
  geom_node_text(aes(label = name),
                 size = node_size_proj/2,
                 repel = TRUE,
                 force = 3,
                 segment.size = 0.2,
                 segment.alpha = 0.5,
                 max.overlaps = Inf,
                 show.legend = FALSE) +
  scale_edge_width_continuous(range = c(0.2, 3)) +
  scale_edge_alpha_continuous(range = c(0.2, 0.7)) +
  scale_color_brewer(palette = "Set1",
                     name = "Comunidad") +
  labs(title = "Red de Proyección de Sustancias",
       subtitle = "Basado en pacientes compartidos | Datos completos") +
  theme_void() +
  theme(plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        plot.subtitle = element_text(size = 12, hjust = 0.5),
        legend.position = "right")

Red de proyección de sustancias (datos completos)

La proyección de la red paciente–sustancia sobre el modo “sustancias” convierte las relaciones indirectas en enlaces directos entre sustancias ponderados por el número de pacientes que comparten, lo que permite analizar la estructura de co‑ocurrencia desde una perspectiva agregada y detectar pares de sustancias que aparecen sistemáticamente en los mismos perfiles de consumo2527. A diferencia de una red de co‑ocurrencia construida solo con episodios simultáneos, esta proyección integra la experiencia compartida de los pacientes a lo largo del período observado, de modo que las asociaciones reflejan patrones de consumo comunes aunque no coincidan en el tiempo, tal como se justifica en el marco de las redes temporales y en aplicaciones de proyección a partir de datos longitudinales28.

4.3 Red Interactiva de Proyección

Código
# Calcular layout con mayor separación
set.seed(123)
coords_proj <- layout_with_fr(g_substances, weights = E(g_substances)$weight)
coords_proj <- coords_proj * 300

# Calcular tamaños proporcionales
node_strength_proj <- strength(g_substances)
node_sizes_proj <- (node_strength_proj / max(node_strength_proj)) * 50 + 10

# Calcular k-core
coreness_bi <- coreness(g_substances)

# Preparar datos para visNetwork
nodes_proj <- data.frame(
  id = V(g_substances)$name,
  label = V(g_substances)$name,
  value = node_sizes_proj,
  group = membership(communities_proj),
  x = coords_proj[,1],
  y = coords_proj[,2],
  title = paste0(
    "<b>", V(g_substances)$name, "</b><br>",
    "Fuerza de conexión: ", format(round(node_strength_proj, 0), big.mark = ","), "<br>",
    "Conexiones: ", degree(g_substances), "<br>",
    "Comunidad: ", membership(communities_proj), "<br>",
    "K-Core: ", coreness_bi
  ),
  font.size = node_sizes_proj * 0.8,  # Tamaño de letra proporcional al nodo
  font.face = "Arial",  # Fuente Arial
  physics = TRUE
)

# Ajustar posiciones de nodos principales
if ("Alcohol" %in% nodes_proj$id) {
  idx <- which(nodes_proj$id == "Alcohol")
  nodes_proj[idx, c("x", "y")] <- c(-400, 0)
  nodes_proj[idx, "fixed"] <- TRUE
}
if ("Marihuana" %in% nodes_proj$id) {
  idx <- which(nodes_proj$id == "Marihuana")
  nodes_proj[idx, c("x", "y")] <- c(400, 0)
  nodes_proj[idx, "fixed"] <- TRUE
}
if ("Pasta Base" %in% nodes_proj$id) {
  idx <- which(nodes_proj$id == "Pasta Base")
  nodes_proj[idx, c("x", "y")] <- c(0, 400)
  nodes_proj[idx, "fixed"] <- TRUE
}
if ("Cocaína" %in% nodes_proj$id) {
  idx <- which(nodes_proj$id == "Cocaína")
  nodes_proj[idx, c("x", "y")] <- c(0, -400)
  nodes_proj[idx, "fixed"] <- TRUE
}

# Preparar edges con grosor proporcional
edges_proj_df <- as_data_frame(g_substances, what = "edges")
edges_proj <- data.frame(
  from = edges_proj_df$from,
  to = edges_proj_df$to,
  value = edges_proj_df$weight / 50,
  title = paste0(
    edges_proj_df$from, " ↔ ", edges_proj_df$to, "<br>",
    "Pacientes compartidos: ", format(edges_proj_df$weight, big.mark = ",")
  ),
  color = "rgba(150,150,150,0.4)",
  highlight = "#FF6B6B"
)

# Crear visualización interactiva
visNetwork(nodes_proj, edges_proj, height = "700px", width = "100%") %>%
  visOptions(
    highlightNearest = list(
      enabled = TRUE,
      hover = TRUE,
      degree = 1,
      labelOnly = FALSE
    ),
    nodesIdSelection = list(
      enabled = TRUE,
      style = 'width: 200px; height: 30px;',
      main = "Buscar sustancia"
    ),
    selectedBy = list(
      variable = "group",
      style = 'width: 200px; height: 30px;',
      main = "Filtrar por comunidad"
    )
  ) %>%
  visPhysics(
    enabled = TRUE,
    stabilization = list(
      enabled = TRUE,
      iterations = 500,
      updateInterval = 50
    ),
    barnesHut = list(
      gravitationalConstant = -10000,
      centralGravity = 0.1,
      springLength = 350,
      springConstant = 0.001,
      damping = 0.5,
      avoidOverlap = 1
    ),
    maxVelocity = 50,
    minVelocity = 0.75,
    solver = "barnesHut"
  ) %>%
  visLayout(
    randomSeed = 123,
    improvedLayout = TRUE
  ) %>%
  visInteraction(
    navigationButtons = TRUE,
    dragNodes = TRUE,
    dragView = TRUE,
    zoomView = TRUE,
    hover = TRUE,
    tooltipDelay = 0,
    hideEdgesOnDrag = TRUE,
    zoomSpeed = 0.5
  ) %>%
  visEdges(
    smooth = list(
      enabled = TRUE,
      type = "continuous",
      roundness = 0.5
    ),
    width = 2,
    physics = TRUE,
    scaling = list(
      min = 0.5,
      max = 6
    ),
    color = list(
      color = "rgba(150,150,150,0.4)",
      highlight = "#FF6B6B",
      hover = "#FF6B6B"
    )
  ) %>%
  visNodes(
    shape = "dot",
    scaling = list(
      min = 10,
      max = 60,
      label = list(
        enabled = TRUE,
        min = 10,
        max = 40,
        drawThreshold = 0,
        maxVisible = 50
      )
    ),
    font = list(
      face = "Arial, sans-serif",
      align = "center"
    ),
    borderWidth = 2,
    borderWidthSelected = 4,
    color = list(
      background = "#e74c3c",
      border = "#c0392b",
      highlight = list(
        background = "#ff6b6b",
        border = "#e74c3c"
      )
    )
  ) %>%
  visLegend(
    enabled = TRUE,
    position = "right",
    width = 0.15,
    useGroups = TRUE,
    main = "Comunidades"
  ) %>%
  visConfigure(enabled = FALSE)

Red interactiva de proyección de sustancias (datos completos)

4.4 Estadísticas Descriptivas de la Red Bipartita Proyectada

Código
# Usar la proyección de la red bipartita como red principal
g_analysis <- g_substances

# Calcular todas las métricas requeridas
metrics_complete <- data.frame(
  Métrica = c("Número de nodos", 
              "Número de enlaces",
              "Densidad",
              "Diámetro", 
              "Longitud de camino medio",
              "Clustering global",
              "Clustering medio",
              "Grado medio",
              "Grado máximo",
              "Grado mínimo",
              "Asortatividad",
              "Modularidad"),
  Valor = c(vcount(g_analysis),
            ecount(g_analysis),
            round(edge_density(g_analysis), 4),
            diameter(g_analysis, weights = NA),
            round(mean_distance(g_analysis, weights = NA), 3),
            round(transitivity(g_analysis, type = "global"), 3),
            round(transitivity(g_analysis, type = "average"), 3),
            round(mean(degree(g_analysis)), 2),
            max(degree(g_analysis)),
            min(degree(g_analysis)),
            round(assortativity_degree(g_analysis), 3),
            round(modularity(communities_proj), 3))
)

metrics_complete %>%
  kable(format = "html", align = c("l", "r")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Estadísticas descriptivas completas de la red bipartita proyectada
Métrica Valor
Número de nodos 22.0000
Número de enlaces 188.0000
Densidad 0.8139
Diámetro 2.0000
Longitud de camino medio 1.1860
Clustering global 0.8840
Clustering medio 0.9070
Grado medio 17.0900
Grado máximo 21.0000
Grado mínimo 6.0000
Asortatividad -0.2600
Modularidad 0.0060

4.5 Medidas de Centralidad Completas (Todas las Sustancias)

Código
# Calcular métricas de centralidad
centrality_metrics_all <- data.frame(
  Sustancia = V(g_analysis)$name,
  Grado = degree(g_analysis),
  Fuerza = round(strength(g_analysis), 0),
  Intermediación = round(betweenness(g_analysis, normalized = TRUE), 4),
  Cercanía = round(closeness(g_analysis, normalized = TRUE), 4),
  Eigenvector = round(eigen_centrality(g_analysis)$vector, 4),
  PageRank = round(page_rank(g_analysis)$vector * 100, 3),  # Multiplicado por 100 para legibilidad
  Comunidad = V(g_analysis)$community
) %>%
  # Agregar información de pacientes totales
  left_join(
    bipartite_edges_full %>%
      group_by(sustancia) %>%
      summarise(Pacientes_Totales = n_distinct(HASH_KEY)),
    by = c("Sustancia" = "sustancia")
  ) %>%
  # Ordenar por Intermediación
  arrange(desc(Intermediación))

# Mostrar todas las sustancias
centrality_metrics_all %>%
  select(Sustancia, Grado, Fuerza, Intermediación, Cercanía, 
         Eigenvector, PageRank, Pacientes_Totales, Comunidad) %>%
  kable(format = "html", row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = FALSE,
                font_size = 12) %>%
  scroll_box(height = "500px")
Medidas de centralidad para todas las sustancias (ordenadas por Intermediación)
Sustancia Grado Fuerza Intermediación Cercanía Eigenvector PageRank Pacientes_Totales Comunidad
Fenilciclidina 8 20 0.4198 0.3684 0.0001 0.691 8 2
Heroína 18 186 0.2770 0.3559 0.0013 0.782 53 2
Hongos 16 374 0.1594 0.3443 0.0027 0.801 126 2
Hipnóticos 16 568 0.1300 0.3134 0.0040 0.801 241 2
Crack 16 610 0.1217 0.3043 0.0048 0.820 193 2
Esteroides Anabólicos 6 28 0.0887 0.2333 0.0002 0.686 10 1
Éxtasis/MDMA 19 1482 0.0818 0.3134 0.0105 1.092 501 2
Pasta Base 21 94691 0.0762 0.2917 0.8131 17.089 48914 1
Metadona 10 52 0.0708 0.2593 0.0003 0.700 19 2
Otros 20 6223 0.0697 0.2561 0.0509 1.858 2703 2
Tranquilizantes 15 435 0.0279 0.2727 0.0032 0.772 184 2
Metanfetaminas y Otros Derivados 18 664 0.0262 0.2958 0.0048 0.827 222 2
Estimulantes 18 1707 0.0131 0.2561 0.0129 1.048 697 2
Anfetaminas 19 2974 0.0024 0.2165 0.0235 1.260 1004 2
Alcohol 21 130582 0.0000 0.1533 1.0000 24.273 87036 1
Marihuana 21 110205 0.0000 0.1243 0.9002 20.592 52768 1
Sedantes 20 13102 0.0000 0.1963 0.1076 3.245 6025 2
Cocaína 21 96619 0.0000 0.1533 0.8179 18.198 47990 1
Opioides 20 2093 0.0000 0.2308 0.0149 1.271 896 2
Inhalables 17 3129 0.0000 0.2692 0.0255 1.236 1065 2
Lsd 18 1297 0.0000 0.2019 0.0095 1.016 420 2
Alucinógenos 18 1101 0.0000 0.2360 0.0082 0.941 372 2

4.6 Medidas de Centralidad Agregadas para Toda la Red

Código
# Calcular estadísticas descriptivas para cada medida de centralidad
centrality_summary <- data.frame(
  Medida = c("Grado", "Fuerza", "Intermediación", "Cercanía", "Eigenvector", "PageRank"),
  Media = c(
    round(mean(centrality_metrics_all$Grado), 2),
    round(mean(centrality_metrics_all$Fuerza), 2),
    round(mean(centrality_metrics_all$Intermediación), 4),
    round(mean(centrality_metrics_all$Cercanía), 4),
    round(mean(centrality_metrics_all$Eigenvector), 4),
    round(mean(centrality_metrics_all$PageRank), 3)
  ),
  Mediana = c(
    round(median(centrality_metrics_all$Grado), 2),
    round(median(centrality_metrics_all$Fuerza), 2),
    round(median(centrality_metrics_all$Intermediación), 4),
    round(median(centrality_metrics_all$Cercanía), 4),
    round(median(centrality_metrics_all$Eigenvector), 4),
    round(median(centrality_metrics_all$PageRank), 3)
  ),
  Desv_Estándar = c(
    round(sd(centrality_metrics_all$Grado), 2),
    round(sd(centrality_metrics_all$Fuerza), 2),
    round(sd(centrality_metrics_all$Intermediación), 4),
    round(sd(centrality_metrics_all$Cercanía), 4),
    round(sd(centrality_metrics_all$Eigenvector), 4),
    round(sd(centrality_metrics_all$PageRank), 3)
  ),
  Mínimo = c(
    round(min(centrality_metrics_all$Grado), 2),
    round(min(centrality_metrics_all$Fuerza), 2),
    round(min(centrality_metrics_all$Intermediación), 4),
    round(min(centrality_metrics_all$Cercanía), 4),
    round(min(centrality_metrics_all$Eigenvector), 4),
    round(min(centrality_metrics_all$PageRank), 3)
  ),
  Máximo = c(
    round(max(centrality_metrics_all$Grado), 2),
    round(max(centrality_metrics_all$Fuerza), 2),
    round(max(centrality_metrics_all$Intermediación), 4),
    round(max(centrality_metrics_all$Cercanía), 4),
    round(max(centrality_metrics_all$Eigenvector), 4),
    round(max(centrality_metrics_all$PageRank), 3)
  ),
  CV = c(
    round(sd(centrality_metrics_all$Grado) / mean(centrality_metrics_all$Grado), 2),
    round(sd(centrality_metrics_all$Fuerza) / mean(centrality_metrics_all$Fuerza), 2),
    round(sd(centrality_metrics_all$Intermediación) / mean(centrality_metrics_all$Intermediación), 2),
    round(sd(centrality_metrics_all$Cercanía) / mean(centrality_metrics_all$Cercanía), 2),
    round(sd(centrality_metrics_all$Eigenvector) / mean(centrality_metrics_all$Eigenvector), 2),
    round(sd(centrality_metrics_all$PageRank) / mean(centrality_metrics_all$PageRank), 2)
  )
)

centrality_summary %>%
  kable(format = "html", 
        col.names = c("Medida", "Media", "Mediana", "Desv. Estándar", 
                     "Mínimo", "Máximo", "Coef. Variación"),
        align = c("l", rep("r", 6))) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Resumen estadístico de las medidas de centralidad para toda la red
Medida Media Mediana Desv. Estándar Mínimo Máximo Coef. Variación
Grado 17.0900 18.0000 4.1500 6.0000 21.0000 0.24
Fuerza 21279.1800 1389.5000 42414.9900 20.0000 130582.0000 1.99
Intermediación 0.0711 0.0271 0.1048 0.0000 0.4198 1.47
Cercanía 0.2566 0.2577 0.0657 0.1243 0.3684 0.26
Eigenvector 0.1735 0.0100 0.3447 0.0001 1.0000 1.99
PageRank 4.5450 1.0320 7.5910 0.6860 24.2730 1.67

4.7 Visualización de la Red según Intermediación

Código
# Visualizar red con tamaño proporcional a intermediación
g_viz <- as_tbl_graph(g_analysis)
betweenness_vals <- betweenness(g_analysis, normalized = TRUE)
node_sizes_betw <- sqrt(betweenness_vals) * 30 + 5

set.seed(42)
ggraph(g_viz, layout = 'fr') +
  geom_edge_link(aes(width = weight, alpha = weight),
                 color = "gray50",
                 show.legend = FALSE) +
  geom_node_point(aes(size = betweenness_vals),
                  color = "#e74c3c",
                  alpha = 0.8) +
  geom_node_text(aes(label = name),
                 size = 3,
                 repel = TRUE,
                 force = 2) +
  scale_size_continuous(range = c(2, 15),
                       name = "Intermediación") +
  scale_edge_width_continuous(range = c(0.2, 2)) +
  scale_edge_alpha_continuous(range = c(0.2, 0.6)) +
  labs(title = "Red con Intermediación como Centralidad",
       subtitle = "Tamaño del nodo proporcional a la intermediación normalizada") +
  theme_void() +
  theme(plot.title = element_text(size = 14, face = "bold"),
        plot.subtitle = element_text(size = 11))

Visualización de la red según intermediación (betweenness)

5 Pregunta Comparativa

¿La red bipartita proyectada de co-ocurrencia de sustancias presenta patrones de formación aleatorios o sigue principios de organización específicos?

5.1 Hipótesis

  1. H0 (Modelo Aleatorio - Erdős-Rényi): La red se forma por asociaciones aleatorias entre sustancias
  2. H1 (Modelo Scale-Free - Barabási-Albert): La red sigue un proceso de adhesión preferencial con distribución de grado de ley de potencias

5.2 Implementación de Modelos Nulos y Bootstrapping

Código
# Parámetros de la red real
n_nodes <- vcount(g_analysis)
n_edges <- ecount(g_analysis)
avg_degree <- mean(degree(g_analysis))
p_connect <- edge_density(g_analysis)

# Número de simulaciones para bootstrapping
n_sims <- 5000

# Función para calcular métricas
calculate_metrics <- function(g) {
  c(
    clustering = transitivity(g, type = "global"),
    path_length = mean_distance(g, directed = FALSE),
    diameter = diameter(g, directed = FALSE),
    assortativity = assortativity_degree(g, directed = FALSE),
    max_degree = max(degree(g)),
    degree_variance = var(degree(g))
  )
}

# Métricas de la red real
real_metrics <- calculate_metrics(g_analysis)

# Simulaciones Bootstrap
set.seed(42)

# 1. Modelo Erdős-Rényi (ER)
er_metrics <- t(replicate(n_sims, {
  g_er <- erdos.renyi.game(n_nodes, p_connect, type = "gnp")
  calculate_metrics(g_er)
}))

# 2. Modelo Barabási-Albert (BA)
ba_metrics <- t(replicate(n_sims, {
  m_ba <- max(1, round(avg_degree/2))
  g_ba <- barabasi.game(n_nodes, m = m_ba, directed = FALSE)
  calculate_metrics(g_ba)
}))

# Calcular intervalos de confianza (95%)
calculate_ci <- function(metrics_matrix, confidence = 0.95) {
  alpha <- 1 - confidence
  lower <- apply(metrics_matrix, 2, quantile, probs = alpha/2, na.rm = TRUE)
  upper <- apply(metrics_matrix, 2, quantile, probs = 1 - alpha/2, na.rm = TRUE)
  mean_val <- colMeans(metrics_matrix, na.rm = TRUE)
  return(list(mean = mean_val, lower = lower, upper = upper))
}

er_ci <- calculate_ci(er_metrics)
ba_ci <- calculate_ci(ba_metrics)

# Crear dataframe con resultados incluyendo intervalos de confianza
results_df <- data.frame(
  Métrica = c("Clustering", "Path Length", "Diameter", 
              "Assortativity", "Max Degree", "Degree Variance"),
  Real = real_metrics,
  ER_mean = er_ci$mean,
  ER_CI_lower = er_ci$lower,
  ER_CI_upper = er_ci$upper,
  BA_mean = ba_ci$mean,
  BA_CI_lower = ba_ci$lower,
  BA_CI_upper = ba_ci$upper
)

results_df %>%
  mutate(across(where(is.numeric), round, 3)) %>%
  mutate(
    ER = paste0(round(ER_mean, 3), " [", round(ER_CI_lower, 3), ", ", round(ER_CI_upper, 3), "]"),
    BA = paste0(round(BA_mean, 3), " [", round(BA_CI_lower, 3), ", ", round(BA_CI_upper, 3), "]")
  ) %>%
  select(Métrica, Real, ER, BA) %>%
  kable(format = "html",
        caption = "Comparación de métricas: Red real vs Modelos nulos (1000 simulaciones, IC 95%)",
        col.names = c("Métrica", "Red Real", "Erdős-Rényi [IC 95%]", "Barabási-Albert [IC 95%]"),
        row.names = F) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Comparación de métricas: Red real vs Modelos nulos (1000 simulaciones, IC 95%)
Métrica Red Real Erdős-Rényi [IC 95%] Barabási-Albert [IC 95%]
Clustering 0.884 0.812 [0.757, 0.862] 0.708 [0.694, 0.722]
Path Length 4.208 1.187 [1.139, 1.238] 1.338 [1.338, 1.338]
Diameter 10.000 2 [2, 2] 2 [2, 2]
Assortativity -0.260 -0.092 [-0.152, -0.029] -0.13 [-0.224, -0.032]
Max Degree 21.000 20.034 [19, 21] 19.383 [18, 21]
Degree Variance 17.229 3.017 [1.515, 5.16] 11.447 [8.656, 14.563]
Código
# Preparar datos para visualización
metrics_long <- data.frame(
  Clustering = c(er_metrics[,1], ba_metrics[,1]),
  PathLength = c(er_metrics[,2], ba_metrics[,2]),
  Assortativity = c(er_metrics[,4], ba_metrics[,4]),
  MaxDegree = c(er_metrics[,5], ba_metrics[,5]),
  Model = rep(c("ER", "BA"), each = n_sims)
)

# Gráficos de distribución con intervalos de confianza como líneas verticales
p1 <- ggplot(metrics_long, aes(x = Clustering, fill = Model)) +
  geom_density(alpha = 0.5) +
  # Línea vertical para el valor real
  geom_vline(xintercept = real_metrics[1], color = "red", linetype = "dashed", size = 1.2) +
  # Intervalos de confianza como líneas verticales para ER
  geom_vline(xintercept = er_ci$lower[1], color = "#3498db", linetype = "dotted", size = 1) +
  geom_vline(xintercept = er_ci$upper[1], color = "#3498db", linetype = "dotted", size = 1) +
  # Intervalos de confianza como líneas verticales para BA
  geom_vline(xintercept = ba_ci$lower[1], color = "#9b59b6", linetype = "dotted", size = 1) +
  geom_vline(xintercept = ba_ci$upper[1], color = "#9b59b6", linetype = "dotted", size = 1) +
  scale_fill_manual(values = c("#3498db", "#9b59b6")) +
  labs(title = "Clustering", x = "Valor", y = "Densidad") +
  theme_minimal() +
  annotate("text", x = real_metrics[1], y = 0, label = "Real", color = "red", vjust = -1)

p2 <- ggplot(metrics_long, aes(x = PathLength, fill = Model)) +
  geom_density(alpha = 0.5) +
  geom_vline(xintercept = real_metrics[2], color = "red", linetype = "dashed", size = 1.2) +
  # Intervalos de confianza como líneas verticales para ER
  geom_vline(xintercept = er_ci$lower[2], color = "#3498db", linetype = "dotted", size = 1) +
  geom_vline(xintercept = er_ci$upper[2], color = "#3498db", linetype = "dotted", size = 1) +
  # Intervalos de confianza como líneas verticales para BA
  geom_vline(xintercept = ba_ci$lower[2], color = "#9b59b6", linetype = "dotted", size = 1) +
  geom_vline(xintercept = ba_ci$upper[2], color = "#9b59b6", linetype = "dotted", size = 1) +
  scale_fill_manual(values = c("#3498db", "#9b59b6")) +
  labs(title = "Path Length", x = "Valor", y = "Densidad") +
  theme_minimal() +
  annotate("text", x = real_metrics[2], y = 0, label = "Real", color = "red", vjust = -1)

p3 <- ggplot(metrics_long, aes(x = Assortativity, fill = Model)) +
  geom_density(alpha = 0.5) +
  geom_vline(xintercept = real_metrics[4], color = "red", linetype = "dashed", size = 1.2) +
  # Intervalos de confianza como líneas verticales para ER
  geom_vline(xintercept = er_ci$lower[4], color = "#3498db", linetype = "dotted", size = 1) +
  geom_vline(xintercept = er_ci$upper[4], color = "#3498db", linetype = "dotted", size = 1) +
  # Intervalos de confianza como líneas verticales para BA
  geom_vline(xintercept = ba_ci$lower[4], color = "#9b59b6", linetype = "dotted", size = 1) +
  geom_vline(xintercept = ba_ci$upper[4], color = "#9b59b6", linetype = "dotted", size = 1) +
  scale_fill_manual(values = c("#3498db", "#9b59b6")) +
  labs(title = "Assortativity", x = "Valor", y = "Densidad") +
  theme_minimal() +
  annotate("text", x = real_metrics[4], y = 0, label = "Real", color = "red", vjust = -1)

p4 <- ggplot(metrics_long, aes(x = MaxDegree, fill = Model)) +
  geom_density(alpha = 0.5) +
  geom_vline(xintercept = real_metrics[5], color = "red", linetype = "dashed", size = 1.2) +
  # Intervalos de confianza como líneas verticales para ER
  geom_vline(xintercept = er_ci$lower[5], color = "#3498db", linetype = "dotted", size = 1) +
  geom_vline(xintercept = er_ci$upper[5], color = "#3498db", linetype = "dotted", size = 1) +
  # Intervalos de confianza como líneas verticales para BA
  geom_vline(xintercept = ba_ci$lower[5], color = "#9b59b6", linetype = "dotted", size = 1) +
  geom_vline(xintercept = ba_ci$upper[5], color = "#9b59b6", linetype = "dotted", size = 1) +
  scale_fill_manual(values = c("#3498db", "#9b59b6")) +
  labs(title = "Max Degree", x = "Valor", y = "Densidad") +
  theme_minimal() +
  annotate("text", x = real_metrics[5], y = 0, label = "Real", color = "red", vjust = -1)

# Combinar los gráficos
(p1 + p2) / (p3 + p4) + 
  plot_layout(guides = "collect") & 
  theme(legend.position = "bottom")

Distribuciones bootstrap de métricas clave con intervalos de confianza
Código
# Generar distribuciones de grado con IC para datos reales
set.seed(42)
n_bootstrap <- 5000

# Función para calcular distribución de grado
calc_degree_dist <- function(g) {
  deg <- degree(g)
  deg_table <- table(deg)
  data.frame(
    k = as.numeric(names(deg_table)),
    p_k = as.numeric(deg_table) / sum(deg_table)
  )
}

# Bootstrap para la red real - resampleando nodos
real_dist_bootstrap <- list()
nodes_total <- vcount(g_analysis)

for(i in 1:n_bootstrap) {
  # Resamplear nodos con reemplazo
  sampled_nodes <- sample(1:nodes_total, nodes_total, replace = TRUE)
  g_sampled <- induced_subgraph(g_analysis, sampled_nodes)
  
  # Calcular distribución del grafo resampleado
  if(vcount(g_sampled) > 0 && ecount(g_sampled) > 0) {
    real_dist_bootstrap[[i]] <- calc_degree_dist(g_sampled)
  }
}

# Calcular distribución promedio e IC para la red real
aggregate_real_distribution <- function(dist_list) {
  all_k <- sort(unique(unlist(lapply(dist_list, function(x) x$k))))
  
  # Crear matriz para almacenar todas las p_k
  n_valid <- length(dist_list)
  p_k_matrix <- matrix(0, nrow = n_valid, ncol = length(all_k))
  colnames(p_k_matrix) <- as.character(all_k)
  
  for(i in 1:n_valid) {
    dist <- dist_list[[i]]
    for(j in 1:nrow(dist)) {
      col_idx <- which(all_k == dist$k[j])
      p_k_matrix[i, col_idx] <- dist$p_k[j]
    }
  }
  
  # Calcular media y desviación estándar
  mean_p_k <- colMeans(p_k_matrix, na.rm = TRUE)
  sd_p_k <- apply(p_k_matrix, 2, sd, na.rm = TRUE)
  
  # IC usando ±1 desviación estándar (como en el ejemplo del paper)
  data.frame(
    k = all_k,
    p_k = mean_p_k,
    p_k_lower = pmax(mean_p_k - sd_p_k, 0),  # No puede ser negativo
    p_k_upper = mean_p_k + sd_p_k,
    Model = "Real"
  )
}

# Calcular IC para datos reales
dd_real_ci <- aggregate_real_distribution(real_dist_bootstrap)

# Generar distribuciones teóricas para ER y BA (sin IC, solo línea teórica)
# Para ER: usar la fórmula teórica de distribución binomial
n_nodes <- vcount(g_analysis)
p_connect <- edge_density(g_analysis)
k_range <- 0:50

# Distribución teórica ER (binomial)
dd_er_theory <- data.frame(
  k = k_range,
  p_k = dbinom(k_range, n_nodes - 1, p_connect),
  p_k_lower = NA,
  p_k_upper = NA,
  Model = "ER"
)

# Para BA: generar una muestra grande para obtener la distribución esperada
m_ba <- max(1, round(mean(degree(g_analysis))/2))
g_ba_large <- barabasi.game(n_nodes * 10, m = m_ba, directed = FALSE)
dd_ba_temp <- calc_degree_dist(g_ba_large)

# Interpolar para tener valores en el mismo rango de k
dd_ba_theory <- data.frame(
  k = k_range,
  p_k = approx(dd_ba_temp$k, dd_ba_temp$p_k, xout = k_range, rule = 2)$y,
  p_k_lower = NA,
  p_k_upper = NA,
  Model = "BA"
)
dd_ba_theory$p_k[is.na(dd_ba_theory$p_k)] <- 0

# Combinar todos los datos
dd_all <- rbind(dd_real_ci, dd_er_theory, dd_ba_theory)

# Escala lineal con IC solo para datos reales
p1 <- ggplot(dd_all, aes(x = k, y = p_k, color = Model)) +
  # Banda de confianza solo para datos reales
  geom_ribbon(data = dd_real_ci,
              aes(ymin = p_k_lower, ymax = p_k_upper, fill = Model), 
              alpha = 0.3, color = NA) +
  # Líneas para todos los modelos
  geom_line(size = 1.2, alpha = 0.9) +
  scale_color_manual(values = c("Real" = "#e74c3c", "ER" = "#2ecc71", "BA" = "#9b59b6"),
                    labels = c("Real" = "Datos reales (media)", 
                              "ER" = "ER (teórico)", 
                              "BA" = "BA (teórico)")) +
  scale_fill_manual(values = c("Real" = "#e74c3c"),
                   labels = c("Real" = "Datos reales (±1 sd)"),
                   guide = guide_legend(override.aes = list(alpha = 0.3))) +
  labs(title = "(a) Escala lineal",
       x = "k", y = expression(p[k])) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold"),
        legend.position = "bottom",
        legend.title = element_blank()) +
  coord_cartesian(xlim = c(0, 100))

# Escala log-log con IC solo para datos reales
# Filtrar valores positivos
dd_log <- dd_all %>%
  filter(k > 0, p_k > 0)

dd_real_log <- dd_real_ci %>%
  filter(k > 0, p_k > 0) %>%
  mutate(
    p_k_lower_safe = pmax(p_k_lower, 1e-6),
    p_k_upper_safe = p_k_upper
  )

p2 <- ggplot(dd_log, aes(x = k, y = p_k, color = Model)) +
  # Banda de confianza solo para datos reales
  geom_ribbon(data = dd_real_log,
              aes(ymin = p_k_lower_safe, ymax = p_k_upper_safe, fill = Model), 
              alpha = 0.3, color = NA) +
  # Líneas para todos los modelos
  geom_line(size = 1.2, alpha = 0.9) +
  # Agregar línea de tendencia para BA (ley de potencia)
  geom_smooth(data = dd_log[dd_log$Model == "BA" & dd_log$k > 10,],
              method = "lm", formula = y ~ x, se = FALSE,
              linetype = "dashed", size = 0.8, alpha = 0.5) +
  scale_x_log10(breaks = c(1, 10, 30, 50),
                labels = c("1", "10", "30", "50")) +
  scale_y_log10(breaks = c(1e-4, 1e-3, 1e-2, 1e-1, 1),
                labels = scales::trans_format("log10", scales::math_format(10^.x))) +
  scale_color_manual(values = c("Real" = "#e74c3c", "ER" = "#2ecc71", "BA" = "#9b59b6"),
                    labels = c("Real" = "Datos reales (media)", 
                              "ER" = "ER (teórico)", 
                              "BA" = "BA (teórico)")) +
  scale_fill_manual(values = c("Real" = "#e74c3c"),
                   labels = c("Real" = "Datos reales (±1 sd)"),
                   guide = guide_legend(override.aes = list(alpha = 0.3))) +
  labs(title = "(b) Escala log-log",
       x = "k (log)", y = expression(p[k]~"(log)")) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold"),
        legend.position = "bottom",
        legend.title = element_blank(),
        panel.grid.minor = element_line(size = 0.1))

# Combinar ambos gráficos
p1 + p2 + 
  plot_layout(guides = "collect") & 
  theme(legend.position = "bottom")

Comparación de distribuciones de grado con modelos teóricos e intervalos de confianza

5.3 Tests Estadísticos

Código
# Función para calcular z-scores
calculate_z_scores <- function(real_val, sim_vals) {
  z <- (real_val - mean(sim_vals, na.rm = TRUE)) / sd(sim_vals, na.rm = TRUE)
  p_val <- 2 * pnorm(-abs(z))
  return(c(z_score = z, p_value = p_val))
}

# Calcular z-scores para cada métrica y modelo
test_results <- data.frame(
  Métrica = c("Agrupamiento", "Longitud de camino", "Diámetro", 
              "Asortatividad", "Grado máximo", "Varianza de grado")
)

for (i in 1:length(real_metrics)) {
  # ER
  er_test <- calculate_z_scores(real_metrics[i], er_metrics[,i])
  test_results[i, "ER_Z"] <- er_test[1]
  test_results[i, "ER_p"] <- er_test[2]
  
  # BA
  ba_test <- calculate_z_scores(real_metrics[i], ba_metrics[,i])
  test_results[i, "BA_Z"] <- ba_test[1]
  test_results[i, "BA_p"] <- ba_test[2]
}

# Agregar significancia
test_results <- test_results %>%
  mutate(
    ER_sig = case_when(
      ER_p < 0.001 ~ "***",
      ER_p < 0.01 ~ "**",
      ER_p < 0.05 ~ "*",
      TRUE ~ "ns"
    ),
    BA_sig = case_when(
      BA_p < 0.001 ~ "***",
      BA_p < 0.01 ~ "**",
      BA_p < 0.05 ~ "*",
      TRUE ~ "ns"
    )
  )

test_results %>%
  mutate(across(where(is.numeric), round, 3)) %>%
  kable(format = "html",
        col.names = c("Métrica", "ER Z", "ER p", "BA Z", "BA p", "ER sig", "BA sig"),
        caption = "Z-scores y valores p: *** p<0.001, ** p<0.01, * p<0.05, ns: no significativo") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Tests estadísticos de comparación con modelos nulos
Métrica ER Z ER p BA Z BA p ER sig BA sig
Agrupamiento 2.700 0.007 23.812 0.000 ** ***
Longitud de camino 115.999 0.000 Inf 0.000 *** ***
Diámetro Inf 0.000 Inf 0.000 *** ***
Asortatividad -5.402 0.000 -2.618 0.009 *** **
Grado máximo 1.381 0.167 2.040 0.041 ns *
Varianza de grado 15.058 0.000 3.763 0.000 *** ***

5.4 Interpretación y Conclusiones

Código
# Crear tabla resumen
summary_df <- data.frame(
  Modelo = c("Erdős-Rényi (ER)", "Barabási-Albert (BA)"),
  `Clustering_Fit` = c(
    ifelse(abs(test_results$ER_Z[1]) < 2, "✓", "✗"),
    ifelse(abs(test_results$BA_Z[1]) < 2, "✓", "✗")
  ),
  `PathLength_Fit` = c(
    ifelse(abs(test_results$ER_Z[2]) < 2, "✓", "✗"),
    ifelse(abs(test_results$BA_Z[2]) < 2, "✓", "✗")
  ),
  `Assortativity_Fit` = c(
    ifelse(abs(test_results$ER_Z[4]) < 2, "✓", "✗"),
    ifelse(abs(test_results$BA_Z[4]) < 2, "✓", "✗")
  ),
  `Degree_Dist_Fit` = c(
    ifelse(abs(test_results$ER_Z[5]) < 2 & abs(test_results$ER_Z[6]) < 2, "✓", "✗"),
    ifelse(abs(test_results$BA_Z[5]) < 2 & abs(test_results$BA_Z[6]) < 2, "✓", "✗")
  ),
  `Ajuste_Global` = c("Bajo", "Moderado-Alto"),
  `Conclusión` = c(
    "No captura estructura",
    "Captura heterogeneidad de grado"
  )
)

summary_df %>%
  kable(format = "html",
        col.names = c("Modelo", "Clustering", "Path Length", 
                     "Asortatividad", "Dist. Grado", 
                     "Ajuste Global", "Conclusión")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Resumen de ajuste de modelos
Modelo Clustering Path Length Asortatividad Dist. Grado Ajuste Global Conclusión
Erdős-Rényi (ER) Bajo No captura estructura
Barabási-Albert (BA) Moderado-Alto Captura heterogeneidad de grado

Referencias

1. Volkow, N. D., & Blanco, C. (2023). Substance use disorders: a comprehensive update of classification, epidemiology, neurobiology, clinical aspects, treatment and prevention. World Psychiatry, 22(2), 203-229. https://doi.org/10.1002/wps.21073
2. Connery, H. S., McHugh, R. K., Reilly, M., Shin, S., & Greenfield, S. F. (2020). Substance Use Disorders in Global Mental Health Delivery: Epidemiology, Treatment Gap, and Implementation of Evidence-Based Treatments. Harvard Review of Psychiatry, 28(5), 316-327. https://doi.org/10.1097/HRP.0000000000000271
3. Gómez-Sánchez-Lafuente, C., Guzman-Parra, J., Suarez-Perez, J., Bordallo-Aragon, A., Rodriguez-de-Fonseca, F., & Mayoral-Cleries, F. (2022). Trends in Psychiatric Hospitalizations of Patients With Dual Diagnosis in Spain. Journal of Dual Diagnosis, 18(2), 92-100. https://doi.org/10.1080/15504263.2022.2053770
4. Saxena, S., Thornicroft, G., Knapp, J., & Whiteford, M. (2007). Resources for mental health: scarcity, inequity, and inefficiency. World Psychiatry, 6(1), 1-10. https://doi.org/10.1002/wps.21073
5. Gómez-Sánchez-Lafuente, C., Guzman-Parra, J., Suarez-Perez, J., Mayoral-Cleries, F., & Rodriguez-de-Fonseca, F. (2016). Dual Diagnosis in Spain: Prevalence, Sociodemographic Profile, and Psychiatric Comorbidity in a Sample of Patients Admitted to Psychiatric Inpatient Units. Journal of Dual Diagnosis, 12(3-4), 249-258. https://doi.org/10.1080/15504263.2016.1220207
6. Rojas, G., Gaete, M., Guajardo, M., Martínez, M., Martínez, M., Fritsch, R., & Araya, R. (2002). Prevalencia de trastornos psiquiátricos en pacientes hospitalizados en un hospital general. Revista Médica de Chile, 130(6), 689-696. https://doi.org/10.4067/S0034-98872002000600008
7. Jones, P. J., Ma, R., & McNally, R. J. (2021). Bridge Centrality: A Network Approach to Understanding Comorbidity. Multivariate Behavioral Research, 56(2), 353-367. https://doi.org/10.1080/00273171.2019.1614898
8. Borsboom, D. (2017). A Network Theory of Mental Disorders. World Psychiatry, 16(1), 5-13. https://doi.org/10.1002/wps.20375
9. López-Toro, E., Wolf, C. J. H., González, R. A., Brink, W. van den, Schellekens, A., & Vélez-Pastrana, M. C. (2022). Network Analysis of DSM Symptoms of Substance Use Disorders and Frequently Co-Occurring Mental Disorders in Patients with Substance Use Disorder Who Seek Treatment. Journal of Clinical Medicine, 11(10), 2883. https://doi.org/10.3390/jcm11102883
10. Sun, E. C., Dixit, A., Humphreys, K., Darnall, B. D., Baker, L. C., & Mackey, S. (2017). Association between concurrent use of prescription opioids and benzodiazepines and overdose: retrospective analysis. BMJ, 356, j760. https://doi.org/10.1136/bmj.j760
11. Hernandez, I., He, M., Brooks, M. M., & Zhang, Y. (2018). Exposure-Response Association Between Concurrent Opioid and Benzodiazepine Use and Risk of Opioid-Related Overdose in Medicare Part D Beneficiaries. JAMA Network Open, 1(2), e180919. https://doi.org/10.1001/jamanetworkopen.2018.0919
12. Liu, S., O’Donnell, J., Gladden, R. M., McGlone, L., & Chowdhury, F. (2021). Trends in Nonfatal and Fatal Overdoses Involving Benzodiazepines — 38 States and the District of Columbia, 2019–2020. MMWR. Morbidity and Mortality Weekly Report, 70(34), 1136-1141. https://doi.org/10.15585/mmwr.mm7034a2
13. Pennings, E. J. M., Leccese, A. P., & Wolff, F. A. de. (2002). Effects of Concurrent Use of Alcohol and Cocaine. Addiction, 97(7), 773-783. https://doi.org/10.1046/j.1360-0443.2002.00158.x
14. Mattson, C. L., Tanz, L. J., Quinn, K., Kariisa, M., Patel, R., & Davis, N. L. (2021). Trends and Geographic Patterns in Drug and Synthetic Opioid Involvement in Overdose Deaths — United States, 2013–2019. MMWR. Morbidity and Mortality Weekly Report, 70(6), 202-207. https://doi.org/10.15585/mmwr.mm7006a4
15. Clauset, A., Shalizi, C. R., & Newman, M. E. J. (2009). Power-Law Distributions in Empirical Data. SIAM Review, 51(4), 661-703. https://doi.org/10.1137/070710111
16. Serrano, M. Á., Boguña, M., & Vespignani, A. (2009). Extracting the Multiscale Backbone of Complex Weighted Networks. Proceedings of the National Academy of Sciences of the United States of America, 106(16), 6483-6488. https://doi.org/10.1073/pnas.0808904106
17. Wang, T., Bendayan, R., Msosa, Y., Pritchard, M., Roberts, A., Stewart, R., & Dobson, R. (2022). Patient-centric characterization of multimorbidity trajectories in patients with severe mental illnesses: A temporal bipartite network modeling approach. Journal of Biomedical Informatics, 127, 104010. https://doi.org/10.1016/j.jbi.2022.104010
18. Borgatti, S. P., & Everett, M. G. (1997). Network analysis of 2-mode data. Social Networks, 19(3), 243-269. https://doi.org/10.1016/S0378-8733(96)00301-2
19. Yang, K.-C., Aronson, B., Odabas, M., Ahn, Y.-Y., & Perry, B. L. (2022). Comparing measures of centrality in bipartite patient–prescriber networks: A study of drug seeking for opioid analgesics. PLOS ONE, 17(8), e0273569. https://doi.org/10.1371/journal.pone.0273569
20. Pavlopoulos, G. A., Kontou, P. I., Pavlopoulou, A., Bouyioukos, C., Markou, E., & Bagos, P. G. (2018). Bipartite graphs in systems biology and medicine: a survey of methods and applications. GigaScience, 7(4), giy014. https://doi.org/10.1093/gigascience/giy014
21. Dormann, C. F., & Strauss, R. (2014). A method for detecting modules in quantitative bipartite networks. Methods in Ecology and Evolution, 5(1), 90-98. https://doi.org/10.1111/2041-210X.12139
22. Beckett, S. J. (2016). Improved community detection in weighted bipartite networks. Royal Society Open Science, 3(1), 140536. https://doi.org/10.1098/rsos.140536
23. Wu, Y., Cao, N., Archambault, D., Shen, Q., Qu, H., & Cui, W. (2017). Evaluation of Graph Sampling: A Visualization Perspective. IEEE Transactions on Visualization and Computer Graphics, 23(1), 401-410. https://doi.org/10.1109/TVCG.2016.2598867
24. Hu, P., & Lau, W. C. (2013). A Survey and Taxonomy of Graph Sampling. arXiv preprint. https://arxiv.org/abs/1308.5865
25. Latapy, M., Magnien, C., & Del Vecchio, N. (2008). Basic notions for the analysis of large two-mode networks. Social Networks, 30(1), 31-48. https://doi.org/10.1016/j.socnet.2007.04.006
26. Newman, M. E. J. (2001). Scientific collaboration networks. I. Network construction and fundamental results. Physical Review E, 64(1), 016131. https://doi.org/10.1103/PhysRevE.64.016131
27. Guillaume, J.-L., & Latapy, M. (2006). Bipartite graphs as models of complex networks. Physica A: Statistical Mechanics and its Applications, 371(2), 795-813. https://doi.org/10.1016/j.physa.2006.04.047
28. Holme, P., & Saramäki, J. (2012). Temporal Networks. Physics Reports, 519(3), 97-125. https://doi.org/10.1016/j.physrep.2012.03.001